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

using namespace std;
namespace OpenMD {

  bool Sticky::initialized_ = false;
  ForceField* Sticky::forceField_ = NULL;
  map<int, AtomType*> Sticky::StickyMap;
  map<pair<AtomType*, AtomType*>, StickyInteractionData> Sticky::MixingMap;
  
  Sticky* Sticky::_instance = NULL;
  
  Sticky* Sticky::Instance() {
    if (!_instance) {
      _instance = new Sticky();
    }
    return _instance;
  }
  
  StickyParam Sticky::getStickyParam(AtomType* atomType) {
    
    // Do sanity checking on the AtomType we were passed before
    // building any data structures:
    if (!atomType->isSticky() && !atomType->isStickyPower()) {
      sprintf( painCave.errMsg,
               "Sticky::getStickyParam was passed an atomType (%s) that does\n"
               "\tnot appear to be a Sticky atom.\n",
               atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();
    }
    
    DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
    GenericData* data = daType->getPropertyByName("Sticky");
    if (data == NULL) {
      sprintf( painCave.errMsg, "Sticky::getStickyParam could not find\n"
               "\tSticky parameters for atomType %s.\n", 
               daType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError(); 
    }
    
    StickyParamGenericData* stickyData = dynamic_cast<StickyParamGenericData*>(data);
    if (stickyData == NULL) {
      sprintf( painCave.errMsg,
               "Sticky::getStickyParam could not convert GenericData to\n"
               "\tStickyParamGenericData for atom type %s\n", 
               daType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();          
    }
    
    return stickyData->getData();
  }
  
  void Sticky::initialize() {    
    
    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
    ForceField::AtomTypeContainer::MapTypeIterator i;
    AtomType* at;

    // Sticky handles all of the Sticky-Sticky interactions

    for (at = atomTypes->beginType(i); at != NULL;
         at = atomTypes->nextType(i)) {
      
      if (at->isSticky() || at->isStickyPower())
        addType(at);
    }
    
    initialized_ = true;
  }
      
  void Sticky::addType(AtomType* atomType){
    // add it to the map:
    AtomTypeProperties atp = atomType->getATP();    
    
    pair<map<int,AtomType*>::iterator,bool> ret;    
    ret = StickyMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
    if (ret.second == false) {
      sprintf( painCave.errMsg,
               "Sticky already had a previous entry with ident %d\n",
               atp.ident);
      painCave.severity = OPENMD_INFO;
      painCave.isFatal = 0;
      simError();         
    }
    
    RealType w0i, v0i, v0pi, rli, rui, rlpi, rupi;
    
    StickyParam sticky1 = getStickyParam(atomType);

    // Now, iterate over all known types and add to the mixing map:
    
    map<int, AtomType*>::iterator it;
    for( it = StickyMap.begin(); it != StickyMap.end(); ++it) {
      
      AtomType* atype2 = (*it).second;
    
      StickyParam sticky2 = getStickyParam(atype2);

      StickyInteractionData mixer;        
           
      // Mixing two different sticky types is silly, but if you want 2
      // sticky types in your simulation, we'll let you do it with the
      // Lorentz- Berthelot mixing rules (which happen to do the right thing 
      // when atomType and atype2 happen to be the same.
      
      mixer.rl   = 0.5 * ( sticky1.rl + sticky2.rl );
      mixer.ru   = 0.5 * ( sticky1.ru + sticky2.ru );
      mixer.rlp  = 0.5 * ( sticky1.rlp + sticky2.rlp );
      mixer.rup  = 0.5 * ( sticky1.rup + sticky2.rup );
      mixer.rbig = max(mixer.ru, mixer.rup);
      mixer.w0  = sqrt( sticky1.w0   * sticky2.w0  );
      mixer.v0  = sqrt( sticky1.v0   * sticky2.v0  );
      mixer.v0p = sqrt( sticky1.v0p  * sticky2.v0p );
      mixer.isPower = atomType->isStickyPower() && atype2->isStickyPower();

      CubicSpline* s = new CubicSpline();
      s->addPoint(mixer.rl, 1.0);
      s->addPoint(mixer.ru, 0.0);
      mixer.s = s;

      CubicSpline* sp = new CubicSpline();
      sp->addPoint(mixer.rlp, 1.0);
      sp->addPoint(mixer.rup, 0.0);
      mixer.sp = sp;

      
      pair<AtomType*, AtomType*> key1, key2;
      key1 = make_pair(atomType, atype2);
      key2 = make_pair(atype2, atomType);
      
      MixingMap[key1] = mixer;
      if (key2 != key1) {
        MixingMap[key2] = mixer;
      }
    }
  }

  RealType Sticky::getStickyCut(int atid) { 
    if (!initialized_) initialize();
    std::map<int, AtomType*> :: const_iterator it;
    it = StickyMap.find(atid);
    if (it == StickyMap.end()) {
      sprintf( painCave.errMsg,
               "Sticky::getStickyCut could not find atid %d in StickyMap\n",
               (atid));
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();          
    }

    AtomType* atype = it->second;
    return MixingMap[make_pair(atype, atype)].rbig;
  }

 
  void Sticky::calcForce(AtomType* at1, AtomType* at2, Vector3d d, 
                         RealType rij, RealType r2, RealType sw, 
                         RealType &vpair, RealType &pot, 
                         RotMat3x3d A1, RotMat3x3d A2, Vector3d &f1,
                         Vector3d &t1, Vector3d &t2) {

    // This routine does only the sticky portion of the SSD potential
    // [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].
    // The Lennard-Jones and dipolar interaction must be handled separately.
    
    // We assume that the rotation matrices have already been calculated
    // and placed in the A array.
    
    if (!initialized_) initialize();
    
    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
    StickyInteractionData mixer = MixingMap[key];

    RealType w0  = mixer.w0;
    RealType v0  = mixer.v0; 
    RealType v0p = mixer.v0p;
    RealType rl  = mixer.rl; 
    RealType ru  = mixer.ru; 
    RealType rlp = mixer.rlp;
    RealType rup = mixer.rup;
    RealType rbig = mixer.rbig;
    bool isPower = mixer.isPower;

    if (rij <= rbig) {

      RealType r3 = r2 * rij;
      RealType r5 = r3 * r2;
           
      RotMat3x3d A1trans = A1.transpose();
      RotMat3x3d A2trans = A2.transpose();

      // rotate the inter-particle separation into the two different
      // body-fixed coordinate systems:

      Vector3d ri = A1 * d;

      // negative sign because this is the vector from j to i:

      Vector3d rj = -A2 * d;

      RealType xi = ri.x();
      RealType yi = ri.y();
      RealType zi = ri.z();

      RealType xj = rj.x();
      RealType yj = rj.y();
      RealType zj = rj.z();

      RealType xi2 = xi * xi;
      RealType yi2 = yi * yi;
      RealType zi2 = zi * zi;

      RealType xj2 = xj * xj;
      RealType yj2 = yj * yj;
      RealType zj2 = zj * zj;     

      // calculate the switching info. from the splines

      RealType s = 0.0;
      RealType dsdr = 0.0;
      RealType sp = 0.0;
      RealType dspdr = 0.0;

      if (rij < ru) {
        if (rij < rl) {
          s = 1.0;
          dsdr = 0.0;
        } else {          
          // we are in the switching region 

          pair<RealType, RealType> res = mixer.s->getValueAndDerivativeAt(rij);
          s = res.first;
          dsdr = res.second;
        }
      }

      if (rij < rup) {
        if (rij < rlp) {
          sp = 1.0;
          dspdr = 0.0;
        } else {
          // we are in the switching region 

          pair<RealType, RealType> res =mixer.sp->getValueAndDerivativeAt(rij);
          sp = res.first;
          dspdr = res.second;
        }
      }

      RealType wi = 2.0*(xi2-yi2)*zi / r3;
      RealType wj = 2.0*(xj2-yj2)*zj / r3;
      RealType w = wi+wj;


      RealType zif = zi/rij - 0.6;
      RealType zis = zi/rij + 0.8;

      RealType zjf = zj/rij - 0.6;
      RealType zjs = zj/rij + 0.8;

      RealType wip = zif*zif*zis*zis - w0;
      RealType wjp = zjf*zjf*zjs*zjs - w0;
      RealType wp = wip + wjp;

      Vector3d dwi(4.0*xi*zi/r3  - 6.0*xi*zi*(xi2-yi2)/r5, 
                   - 4.0*yi*zi/r3  - 6.0*yi*zi*(xi2-yi2)/r5,
                   2.0*(xi2-yi2)/r3  - 6.0*zi2*(xi2-yi2)/r5);
      
      Vector3d dwj(4.0*xj*zj/r3  - 6.0*xj*zj*(xj2-yj2)/r5,
                   - 4.0*yj*zj/r3  - 6.0*yj*zj*(xj2-yj2)/r5,
                   2.0*(xj2-yj2)/r3  - 6.0*zj2*(xj2-yj2)/r5);

      RealType uglyi = zif*zif*zis + zif*zis*zis;
      RealType uglyj = zjf*zjf*zjs + zjf*zjs*zjs;

      Vector3d dwip(-2.0*xi*zi*uglyi/r3, 
                    -2.0*yi*zi*uglyi/r3,
                    2.0*(1.0/rij - zi2/r3)*uglyi);

      Vector3d dwjp(-2.0*xj*zj*uglyj/r3,
                    -2.0*yj*zj*uglyj/r3,
                    2.0*(1.0/rij - zj2/r3)*uglyj);

      Vector3d dwidu(4.0*(yi*zi2 + 0.5*yi*(xi2-yi2))/r3,
                     4.0*(xi*zi2 - 0.5*xi*(xi2-yi2))/r3,
                     - 8.0*xi*yi*zi/r3);

      Vector3d dwjdu(4.0*(yj*zj2 + 0.5*yj*(xj2-yj2))/r3,
                     4.0*(xj*zj2 - 0.5*xj*(xj2-yj2))/r3,
                     - 8.0*xj*yj*zj/r3);
      
      Vector3d dwipdu(2.0*yi*uglyi/rij,
                      -2.0*xi*uglyi/rij,
                      0.0);

      Vector3d dwjpdu(2.0*yj*uglyj/rij,
                      -2.0*xj*uglyj/rij,
                      0.0);

      if (isPower) {
        RealType frac1 = 0.25;
        RealType frac2 = 0.75;      
        RealType wi2 = wi*wi;
        RealType wj2 = wj*wj;
        // sticky power has no w' function:
        w = frac1 * wi * wi2 + frac2*wi + frac1*wj*wj2 + frac2*wj + v0p; 
        wp = 0.0;
        dwi = frac1*3.0*wi2*dwi + frac2*dwi;
        dwj = frac1*3.0*wj2*dwi + frac2*dwi;
        dwip = V3Zero;
        dwjp = V3Zero;
        dwidu = frac1*3.0*wi2*dwidu + frac2*dwidu;
        dwidu = frac1*3.0*wj2*dwjdu + frac2*dwjdu;
        dwipdu = V3Zero;
        dwjpdu = V3Zero;
        sp = 0.0;
        dspdr = 0.0;
      }

      vpair += 0.5*(v0*s*w + v0p*sp*wp);
      pot += 0.5*(v0*s*w + v0p*sp*wp)*sw;

      // do the torques first since they are easy:
      // remember that these are still in the body-fixed axes

      Vector3d ti = 0.5*sw*(v0*s*dwidu + v0p*sp*dwipdu);
      Vector3d tj = 0.5*sw*(v0*s*dwjdu + v0p*sp*dwjpdu);

      // go back to lab frame using transpose of rotation matrix:

      t1 += A1trans * ti;
      t2 += A2trans * tj;

      // Now, on to the forces:

      // first rotate the i terms back into the lab frame:

      Vector3d radcomi = (v0 * s * dwi + v0p * sp * dwip) * sw;
      Vector3d radcomj = (v0 * s * dwj + v0p * sp * dwjp) * sw;

      Vector3d fii = A1trans * radcomi;
      Vector3d fjj = A2trans * radcomj;
      
      // now assemble these with the radial-only terms:
      
      f1 += 0.5 * ((v0*dsdr*w + v0p*dspdr*wp) * d / rij + fii - fjj);

    }
      
    return;
    
  }

  void Sticky::do_sticky_pair(int *atid1, int *atid2, RealType *d, 
                              RealType *r, RealType *r2, RealType *sw, 
                              RealType *vpair, RealType *pot, RealType *A1, 
                              RealType *A2, RealType *f1, 
                              RealType *t1, RealType *t2) {
    
    if (!initialized_) initialize();
    
    AtomType* atype1 = StickyMap[*atid1];
    AtomType* atype2 = StickyMap[*atid2];
    
    Vector3d disp(d);
    Vector3d frc(f1);
    Vector3d trq1(t1);
    Vector3d trq2(t2);
    RotMat3x3d Ai(A1);
    RotMat3x3d Aj(A2);
   
    calcForce(atype1, atype2, disp, *r, *r2, *sw, *vpair, *pot, 
              Ai, Aj, frc, trq1, trq2);
      
    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 fortranGetStickyCut FC_FUNC(getstickycut, GETSTICKYCUT)
#define fortranDoStickyPair FC_FUNC(do_sticky_pair, DO_STICKY_PAIR)
  
  RealType fortranGetStickyCut(int* atid) {
    return OpenMD::Sticky::Instance()->getStickyCut(*atid);
  }

  void fortranDoStickyPair(int *atid1, int *atid2, RealType *d, RealType *r, 
                       RealType *r2, RealType *sw, RealType *vpair, RealType *pot, 
                       RealType *A1, RealType *A2, RealType *f1, 
                       RealType *t1, RealType *t2){
    
    return OpenMD::Sticky::Instance()->do_sticky_pair(atid1, atid2, d, r, r2, 
                                                      sw, vpair, pot, A1, A2, 
                                                      f1, t1, t2);
  }
}
Revision:	1485
Committed:	Wed Jul 28 19:52:00 2010 UTC (14 years, 9 months ago) by gezelter
File size:	14973 byte(s)
Log Message:	Converting Sticky over to C++
#	User	Rev	Content
1	gezelter	1485	/*
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/Sticky.hpp"
47			#include "nonbonded/LJ.hpp"
48			#include "utils/simError.h"
49
50			using namespace std;
51			namespace OpenMD {
52
53			bool Sticky::initialized_ = false;
54			ForceField* Sticky::forceField_ = NULL;
55			map<int, AtomType*> Sticky::StickyMap;
56			map<pair<AtomType, AtomType>, StickyInteractionData> Sticky::MixingMap;
57
58			Sticky* Sticky::_instance = NULL;
59
60			Sticky* Sticky::Instance() {
61			if (!_instance) {
62			_instance = new Sticky();
63			}
64			return _instance;
65			}
66
67			StickyParam Sticky::getStickyParam(AtomType* atomType) {
68
69			// Do sanity checking on the AtomType we were passed before
70			// building any data structures:
71			if (!atomType->isSticky() && !atomType->isStickyPower()) {
72			sprintf( painCave.errMsg,
73			"Sticky::getStickyParam was passed an atomType (%s) that does\n"
74			"\tnot appear to be a Sticky atom.\n",
75			atomType->getName().c_str());
76			painCave.severity = OPENMD_ERROR;
77			painCave.isFatal = 1;
78			simError();
79			}
80
81			DirectionalAtomType* daType = dynamic_cast<DirectionalAtomType*>(atomType);
82			GenericData* data = daType->getPropertyByName("Sticky");
83			if (data == NULL) {
84			sprintf( painCave.errMsg, "Sticky::getStickyParam could not find\n"
85			"\tSticky parameters for atomType %s.\n",
86			daType->getName().c_str());
87			painCave.severity = OPENMD_ERROR;
88			painCave.isFatal = 1;
89			simError();
90			}
91
92			StickyParamGenericData* stickyData = dynamic_cast<StickyParamGenericData*>(data);
93			if (stickyData == NULL) {
94			sprintf( painCave.errMsg,
95			"Sticky::getStickyParam could not convert GenericData to\n"
96			"\tStickyParamGenericData for atom type %s\n",
97			daType->getName().c_str());
98			painCave.severity = OPENMD_ERROR;
99			painCave.isFatal = 1;
100			simError();
101			}
102
103			return stickyData->getData();
104			}
105
106			void Sticky::initialize() {
107
108			ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
109			ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
110			ForceField::AtomTypeContainer::MapTypeIterator i;
111			AtomType* at;
112
113			// Sticky handles all of the Sticky-Sticky interactions
114
115			for (at = atomTypes->beginType(i); at != NULL;
116			at = atomTypes->nextType(i)) {
117
118			if (at->isSticky() \|\| at->isStickyPower())
119			addType(at);
120			}
121
122			initialized_ = true;
123			}
124
125			void Sticky::addType(AtomType* atomType){
126			// add it to the map:
127			AtomTypeProperties atp = atomType->getATP();
128
129			pair<map<int,AtomType*>::iterator,bool> ret;
130			ret = StickyMap.insert( pair<int, AtomType*>(atp.ident, atomType) );
131			if (ret.second == false) {
132			sprintf( painCave.errMsg,
133			"Sticky already had a previous entry with ident %d\n",
134			atp.ident);
135			painCave.severity = OPENMD_INFO;
136			painCave.isFatal = 0;
137			simError();
138			}
139
140			RealType w0i, v0i, v0pi, rli, rui, rlpi, rupi;
141
142			StickyParam sticky1 = getStickyParam(atomType);
143
144			// Now, iterate over all known types and add to the mixing map:
145
146			map<int, AtomType*>::iterator it;
147			for( it = StickyMap.begin(); it != StickyMap.end(); ++it) {
148
149			AtomType* atype2 = (*it).second;
150
151			StickyParam sticky2 = getStickyParam(atype2);
152
153			StickyInteractionData mixer;
154
155			// Mixing two different sticky types is silly, but if you want 2
156			// sticky types in your simulation, we'll let you do it with the
157			// Lorentz- Berthelot mixing rules (which happen to do the right thing
158			// when atomType and atype2 happen to be the same.
159
160			mixer.rl = 0.5 * ( sticky1.rl + sticky2.rl );
161			mixer.ru = 0.5 * ( sticky1.ru + sticky2.ru );
162			mixer.rlp = 0.5 * ( sticky1.rlp + sticky2.rlp );
163			mixer.rup = 0.5 * ( sticky1.rup + sticky2.rup );
164			mixer.rbig = max(mixer.ru, mixer.rup);
165			mixer.w0 = sqrt( sticky1.w0 * sticky2.w0 );
166			mixer.v0 = sqrt( sticky1.v0 * sticky2.v0 );
167			mixer.v0p = sqrt( sticky1.v0p * sticky2.v0p );
168			mixer.isPower = atomType->isStickyPower() && atype2->isStickyPower();
169
170			CubicSpline* s = new CubicSpline();
171			s->addPoint(mixer.rl, 1.0);
172			s->addPoint(mixer.ru, 0.0);
173			mixer.s = s;
174
175			CubicSpline* sp = new CubicSpline();
176			sp->addPoint(mixer.rlp, 1.0);
177			sp->addPoint(mixer.rup, 0.0);
178			mixer.sp = sp;
179
180
181			pair<AtomType, AtomType> key1, key2;
182			key1 = make_pair(atomType, atype2);
183			key2 = make_pair(atype2, atomType);
184
185			MixingMap[key1] = mixer;
186			if (key2 != key1) {
187			MixingMap[key2] = mixer;
188			}
189			}
190			}
191
192			RealType Sticky::getStickyCut(int atid) {
193			if (!initialized_) initialize();
194			std::map<int, AtomType*> :: const_iterator it;
195			it = StickyMap.find(atid);
196			if (it == StickyMap.end()) {
197			sprintf( painCave.errMsg,
198			"Sticky::getStickyCut could not find atid %d in StickyMap\n",
199			(atid));
200			painCave.severity = OPENMD_ERROR;
201			painCave.isFatal = 1;
202			simError();
203			}
204
205			AtomType* atype = it->second;
206			return MixingMap[make_pair(atype, atype)].rbig;
207			}
208
209
210			void Sticky::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
211			RealType rij, RealType r2, RealType sw,
212			RealType &vpair, RealType &pot,
213			RotMat3x3d A1, RotMat3x3d A2, Vector3d &f1,
214			Vector3d &t1, Vector3d &t2) {
215
216			// This routine does only the sticky portion of the SSD potential
217			// [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].
218			// The Lennard-Jones and dipolar interaction must be handled separately.
219
220			// We assume that the rotation matrices have already been calculated
221			// and placed in the A array.
222
223			if (!initialized_) initialize();
224
225			pair<AtomType, AtomType> key = make_pair(at1, at2);
226			StickyInteractionData mixer = MixingMap[key];
227
228			RealType w0 = mixer.w0;
229			RealType v0 = mixer.v0;
230			RealType v0p = mixer.v0p;
231			RealType rl = mixer.rl;
232			RealType ru = mixer.ru;
233			RealType rlp = mixer.rlp;
234			RealType rup = mixer.rup;
235			RealType rbig = mixer.rbig;
236			bool isPower = mixer.isPower;
237
238			if (rij <= rbig) {
239
240			RealType r3 = r2 * rij;
241			RealType r5 = r3 * r2;
242
243			RotMat3x3d A1trans = A1.transpose();
244			RotMat3x3d A2trans = A2.transpose();
245
246			// rotate the inter-particle separation into the two different
247			// body-fixed coordinate systems:
248
249			Vector3d ri = A1 * d;
250
251			// negative sign because this is the vector from j to i:
252
253			Vector3d rj = -A2 * d;
254
255			RealType xi = ri.x();
256			RealType yi = ri.y();
257			RealType zi = ri.z();
258
259			RealType xj = rj.x();
260			RealType yj = rj.y();
261			RealType zj = rj.z();
262
263			RealType xi2 = xi * xi;
264			RealType yi2 = yi * yi;
265			RealType zi2 = zi * zi;
266
267			RealType xj2 = xj * xj;
268			RealType yj2 = yj * yj;
269			RealType zj2 = zj * zj;
270
271			// calculate the switching info. from the splines
272
273			RealType s = 0.0;
274			RealType dsdr = 0.0;
275			RealType sp = 0.0;
276			RealType dspdr = 0.0;
277
278			if (rij < ru) {
279			if (rij < rl) {
280			s = 1.0;
281			dsdr = 0.0;
282			} else {
283			// we are in the switching region
284
285			pair<RealType, RealType> res = mixer.s->getValueAndDerivativeAt(rij);
286			s = res.first;
287			dsdr = res.second;
288			}
289			}
290
291			if (rij < rup) {
292			if (rij < rlp) {
293			sp = 1.0;
294			dspdr = 0.0;
295			} else {
296			// we are in the switching region
297
298			pair<RealType, RealType> res =mixer.sp->getValueAndDerivativeAt(rij);
299			sp = res.first;
300			dspdr = res.second;
301			}
302			}
303
304			RealType wi = 2.0(xi2-yi2)zi / r3;
305			RealType wj = 2.0(xj2-yj2)zj / r3;
306			RealType w = wi+wj;
307
308
309			RealType zif = zi/rij - 0.6;
310			RealType zis = zi/rij + 0.8;
311
312			RealType zjf = zj/rij - 0.6;
313			RealType zjs = zj/rij + 0.8;
314
315			RealType wip = zifzifzis*zis - w0;
316			RealType wjp = zjfzjfzjs*zjs - w0;
317			RealType wp = wip + wjp;
318
319			Vector3d dwi(4.0xizi/r3 - 6.0xizi*(xi2-yi2)/r5,
320			- 4.0yizi/r3 - 6.0yizi*(xi2-yi2)/r5,
321			2.0(xi2-yi2)/r3 - 6.0zi2*(xi2-yi2)/r5);
322
323			Vector3d dwj(4.0xjzj/r3 - 6.0xjzj*(xj2-yj2)/r5,
324			- 4.0yjzj/r3 - 6.0yjzj*(xj2-yj2)/r5,
325			2.0(xj2-yj2)/r3 - 6.0zj2*(xj2-yj2)/r5);
326
327			RealType uglyi = zifzifzis + zifziszis;
328			RealType uglyj = zjfzjfzjs + zjfzjszjs;
329
330			Vector3d dwip(-2.0xizi*uglyi/r3,
331			-2.0yizi*uglyi/r3,
332			2.0(1.0/rij - zi2/r3)uglyi);
333
334			Vector3d dwjp(-2.0xjzj*uglyj/r3,
335			-2.0yjzj*uglyj/r3,
336			2.0(1.0/rij - zj2/r3)uglyj);
337
338			Vector3d dwidu(4.0(yizi2 + 0.5yi(xi2-yi2))/r3,
339			4.0(xizi2 - 0.5xi(xi2-yi2))/r3,
340			- 8.0xiyi*zi/r3);
341
342			Vector3d dwjdu(4.0(yjzj2 + 0.5yj(xj2-yj2))/r3,
343			4.0(xjzj2 - 0.5xj(xj2-yj2))/r3,
344			- 8.0xjyj*zj/r3);
345
346			Vector3d dwipdu(2.0yiuglyi/rij,
347			-2.0xiuglyi/rij,
348			0.0);
349
350			Vector3d dwjpdu(2.0yjuglyj/rij,
351			-2.0xjuglyj/rij,
352			0.0);
353
354			if (isPower) {
355			RealType frac1 = 0.25;
356			RealType frac2 = 0.75;
357			RealType wi2 = wi*wi;
358			RealType wj2 = wj*wj;
359			// sticky power has no w' function:
360			w = frac1 * wi * wi2 + frac2wi + frac1wjwj2 + frac2wj + v0p;
361			wp = 0.0;
362			dwi = frac13.0wi2dwi + frac2dwi;
363			dwj = frac13.0wj2dwi + frac2dwi;
364			dwip = V3Zero;
365			dwjp = V3Zero;
366			dwidu = frac13.0wi2dwidu + frac2dwidu;
367			dwidu = frac13.0wj2dwjdu + frac2dwjdu;
368			dwipdu = V3Zero;
369			dwjpdu = V3Zero;
370			sp = 0.0;
371			dspdr = 0.0;
372			}
373
374			vpair += 0.5(v0sw + v0psp*wp);
375			pot += 0.5(v0sw + v0pspwp)sw;
376
377			// do the torques first since they are easy:
378			// remember that these are still in the body-fixed axes
379
380			Vector3d ti = 0.5sw(v0sdwidu + v0pspdwipdu);
381			Vector3d tj = 0.5sw(v0sdwjdu + v0pspdwjpdu);
382
383			// go back to lab frame using transpose of rotation matrix:
384
385			t1 += A1trans * ti;
386			t2 += A2trans * tj;
387
388			// Now, on to the forces:
389
390			// first rotate the i terms back into the lab frame:
391
392			Vector3d radcomi = (v0 * s * dwi + v0p * sp * dwip) * sw;
393			Vector3d radcomj = (v0 * s * dwj + v0p * sp * dwjp) * sw;
394
395			Vector3d fii = A1trans * radcomi;
396			Vector3d fjj = A2trans * radcomj;
397
398			// now assemble these with the radial-only terms:
399
400			f1 += 0.5 * ((v0dsdrw + v0pdspdrwp) * d / rij + fii - fjj);
401
402			}
403
404			return;
405
406			}
407
408			void Sticky::do_sticky_pair(int atid1, int atid2, RealType *d,
409			RealType r, RealType r2, RealType *sw,
410			RealType vpair, RealType pot, RealType *A1,
411			RealType A2, RealType f1,
412			RealType t1, RealType t2) {
413
414			if (!initialized_) initialize();
415
416			AtomType* atype1 = StickyMap[*atid1];
417			AtomType* atype2 = StickyMap[*atid2];
418
419			Vector3d disp(d);
420			Vector3d frc(f1);
421			Vector3d trq1(t1);
422			Vector3d trq2(t2);
423			RotMat3x3d Ai(A1);
424			RotMat3x3d Aj(A2);
425
426			calcForce(atype1, atype2, disp, r, r2, sw, vpair, *pot,
427			Ai, Aj, frc, trq1, trq2);
428
429			f1[0] = frc.x();
430			f1[1] = frc.y();
431			f1[2] = frc.z();
432
433			t1[0] = trq1.x();
434			t1[1] = trq1.y();
435			t1[2] = trq1.z();
436
437			t2[0] = trq2.x();
438			t2[1] = trq2.y();
439			t2[2] = trq2.z();
440
441			return;
442			}
443			}
444
445			extern "C" {
446
447			#define fortranGetStickyCut FC_FUNC(getstickycut, GETSTICKYCUT)
448			#define fortranDoStickyPair FC_FUNC(do_sticky_pair, DO_STICKY_PAIR)
449
450			RealType fortranGetStickyCut(int* atid) {
451			return OpenMD::Sticky::Instance()->getStickyCut(*atid);
452			}
453
454			void fortranDoStickyPair(int atid1, int atid2, RealType d, RealType r,
455			RealType r2, RealType sw, RealType vpair, RealType pot,
456			RealType A1, RealType A2, RealType *f1,
457			RealType t1, RealType t2){
458
459			return OpenMD::Sticky::Instance()->do_sticky_pair(atid1, atid2, d, r, r2,
460			sw, vpair, pot, A1, A2,
461			f1, t1, t2);
462			}
463			}