src/nonbonded/Morse.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/Morse.hpp"
#include "utils/simError.h"
#include "types/NonBondedInteractionType.hpp"

using namespace std;

namespace OpenMD {

  bool Morse::initialized_ = false;
  bool Morse::shiftedPot_ = false;
  bool Morse::shiftedFrc_ = false;
  ForceField* Morse::forceField_ = NULL;
  map<int, AtomType*> Morse::MorseMap;
  map<pair<AtomType*, AtomType*>, MorseInteractionData> Morse::MixingMap;
  map<string, MorseInteractionType> Morse::stringToEnumMap_;
  
  Morse* Morse::_instance = NULL;

  Morse* Morse::Instance() {
    if (!_instance) {
      _instance = new Morse();
    }
    return _instance;
  }

  void Morse::initialize() {    

    stringToEnumMap_["shiftedMorse"] = shiftedMorse;
    stringToEnumMap_["repulsiveMorse"] = repulsiveMorse;

    ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
    ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
    NonBondedInteractionType* nbt;

    for (nbt = nbiTypes->beginType(j); nbt != NULL; 
         nbt = nbiTypes->nextType(j)) {
      
      if (nbt->isMorse()) {
        
        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
        
        GenericData* data = nbt->getPropertyByName("Morse");
        if (data == NULL) {
          sprintf( painCave.errMsg, "Morse::initialize could not find\n"
                   "\tMorse parameters for %s - %s interaction.\n", 
                   atypes.first->getName().c_str(),
                   atypes.second->getName().c_str());
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError(); 
        }
        
        MorseData* morseData = dynamic_cast<MorseData*>(data);
        if (morseData == NULL) {
          sprintf( painCave.errMsg,
                   "Morse::initialize could not convert GenericData to\n"
                   "\tMorseData for %s - %s interaction.\n", 
                   atypes.first->getName().c_str(),
                   atypes.second->getName().c_str());
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError();          
        }
        
        MorseParam morseParam = morseData->getData();

        RealType De = morseParam.De;
        RealType Re = morseParam.Re;
        RealType beta = morseParam.beta;
        string interactionType = morseParam.interactionType;

        toUpper(interactionType);
        map<string, MorseInteractionType>::iterator i;
        i = stringToEnumMap_.find(interactionType);
        if (i != stringToEnumMap_.end()) { 
          addExplicitInteraction(atypes.first, atypes.second, 
                                 De, Re, beta, i->second );
        } else {
          sprintf( painCave.errMsg,
                   "Morse::initialize found unknown Morse interaction type\n"
                   "\t(%s) for %s - %s interaction.\n", 
                   morseParam.interactionType.c_str(),
                   atypes.first->getName().c_str(),
                   atypes.second->getName().c_str());
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError();          
        }
      }
    }  
    initialized_ = true;
  }
      
  void Morse::addExplicitInteraction(AtomType* atype1, AtomType* atype2, 
                                     RealType De, RealType Re, RealType beta, 
                                     MorseInteractionType mit) {

    AtomTypeProperties atp = atype1->getATP();    
    MorseMap.insert( pair<int, AtomType*>(atp.ident, atype1) );

    atp = atype2->getATP();    
    MorseMap.insert( pair<int, AtomType*>(atp.ident, atype2) );

    MorseInteractionData mixer;
    mixer.De = De;
    mixer.Re = Re;
    mixer.beta = beta;
    mixer.interactionType = mit;

    pair<AtomType*, AtomType*> key1, key2;
    key1 = make_pair(atype1, atype2);
    key2 = make_pair(atype2, atype1);
    
    MixingMap[key1] = mixer;
    if (key2 != key1) {
      MixingMap[key2] = mixer;
    }    
  }
  
  void Morse::calcForce(AtomType* at1, AtomType* at2, Vector3d d, 
                        RealType r, RealType r2, RealType rcut, RealType sw, 
                        RealType vdwMult, RealType &vpair, RealType &pot, 
                        Vector3d &f1) {

    if (!initialized_) initialize();
    
    RealType myPot = 0.0;
    RealType myPotC = 0.0;
    RealType myDeriv = 0.0;
    RealType myDerivC = 0.0;

    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
    MorseInteractionData mixer = MixingMap[key];

    RealType De = mixer.De;
    RealType Re = mixer.Re;
    RealType beta = mixer.beta;
    MorseInteractionType interactionType = mixer.interactionType;
   
    // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
    
    RealType expt     = -beta*(r - Re);
    RealType expfnc   = exp(expt);
    RealType expfnc2  = expfnc*expfnc;

    RealType exptC = 0.0;
    RealType expfncC = 0.0;
    RealType expfnc2C = 0.0;

    if (Morse::shiftedPot_ || Morse::shiftedFrc_) {
      exptC     = -beta*(rcut - Re);
      expfncC   = exp(exptC);
      expfnc2C  = expfncC*expfncC;
    }


    switch(interactionType) {
    case shiftedMorse : {

      myPot  = De * (expfnc2  - 2.0 * expfnc);
      myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);

      if (Morse::shiftedPot_) {
        myPotC = De * (expfnc2C - 2.0 * expfncC);
        myDerivC = 0.0;
      } else if (Morse::shiftedFrc_) {
        myPotC = De * (expfnc2C - 2.0 * expfncC);
        myDerivC  = 2.0 * De * beta * (expfnc2C - expfnc2C);
        myPotC += myDerivC * (r - rcut);
      } else {
        myPotC = 0.0;
        myDerivC = 0.0;
      }

      break;
    }
    case repulsiveMorse : {

      myPot  = De * expfnc2;
      myDeriv  = -2.0 * De * beta * expfnc2;

      if (Morse::shiftedPot_) {
        myPotC = De * expfnc2C;
        myDerivC = 0.0;
      } else if (Morse::shiftedFrc_) {
        myPotC = De * expfnc2C;
        myDerivC = -2.0 * De * beta * expfnc2C;
        myPotC += myDerivC * (r - rcut);
      } else {
        myPotC = 0.0;
        myDerivC = 0.0;
      }

      break;
    }
    }

    RealType pot_temp = vdwMult * (myPot - myPotC);
    vpair += pot_temp;

    RealType dudr = sw * vdwMult * (myDeriv - myDerivC);
    
    pot += sw * pot_temp;
    f1 = d * dudr / r;

    return;

  }

  void Morse::do_morse_pair(int *atid1, int *atid2, RealType *d, 
                            RealType *rij, RealType *r2, RealType *rcut, 
                            RealType *sw, RealType *vdwMult,
                            RealType *vpair, RealType *pot, RealType *f1) {

    if (!initialized_) initialize();
    
    AtomType* atype1 = MorseMap[*atid1];
    AtomType* atype2 = MorseMap[*atid2];
    
    Vector3d disp(d[0], d[1], d[2]);
    Vector3d frc(f1[0], f1[1], f1[2]);
    
    calcForce(atype1, atype2, disp, *rij, *r2, *rcut, *sw, *vdwMult, *vpair, 
              *pot, frc);
      
    f1[0] = frc.x();
    f1[1] = frc.y();
    f1[2] = frc.z();

    return;    
  }
    
  void Morse::setMorseDefaultCutoff(RealType *thisRcut, int *sP, int *sF) {
    shiftedPot_ = (bool)(*sP);
    shiftedFrc_ = (bool)(*sF);
  }
}

extern "C" {
  
#define fortranSetMorseCutoff FC_FUNC(setmorsedefaultcutoff, SETMORSEDEFAULTCUTOFF)
#define fortranDoMorsePair FC_FUNC(do_morse_pair, DO_MORSE_PAIR)
  
  void fortranSetMorseCutoff(RealType *rcut, int *shiftedPot, int *shiftedFrc) {
    return OpenMD::Morse::Instance()->setMorseDefaultCutoff(rcut, shiftedPot, 
                                                            shiftedFrc);
  }
  void fortranDoMorsePair(int *atid1, int *atid2, RealType *d, RealType *rij, 
                          RealType *r2, RealType *rcut, RealType *sw, 
                          RealType *vdwMult, RealType* vpair, RealType* pot, 
                          RealType *f1){
    
    return OpenMD::Morse::Instance()->do_morse_pair(atid1, atid2, d, rij, r2, 
                                                    rcut, sw, vdwMult, vpair, 
                                                    pot, f1);
  }
}
Revision:	1501
Committed:	Wed Sep 15 19:32:10 2010 UTC (14 years, 10 months ago) by gezelter
File size:	10196 byte(s)
Log Message:	Starting migration of Morse to C++
#	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/Morse.hpp"
47			#include "utils/simError.h"
48			#include "types/NonBondedInteractionType.hpp"
49
50			using namespace std;
51
52			namespace OpenMD {
53
54			bool Morse::initialized_ = false;
55			bool Morse::shiftedPot_ = false;
56			bool Morse::shiftedFrc_ = false;
57			ForceField* Morse::forceField_ = NULL;
58			map<int, AtomType*> Morse::MorseMap;
59			map<pair<AtomType, AtomType>, MorseInteractionData> Morse::MixingMap;
60			map<string, MorseInteractionType> Morse::stringToEnumMap_;
61
62			Morse* Morse::_instance = NULL;
63
64			Morse* Morse::Instance() {
65			if (!_instance) {
66			_instance = new Morse();
67			}
68			return _instance;
69			}
70
71			void Morse::initialize() {
72
73			stringToEnumMap_["shiftedMorse"] = shiftedMorse;
74			stringToEnumMap_["repulsiveMorse"] = repulsiveMorse;
75
76			ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
77			ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
78			NonBondedInteractionType* nbt;
79
80			for (nbt = nbiTypes->beginType(j); nbt != NULL;
81			nbt = nbiTypes->nextType(j)) {
82
83			if (nbt->isMorse()) {
84
85			pair<AtomType, AtomType> atypes = nbt->getAtomTypes();
86
87			GenericData* data = nbt->getPropertyByName("Morse");
88			if (data == NULL) {
89			sprintf( painCave.errMsg, "Morse::initialize could not find\n"
90			"\tMorse parameters for %s - %s interaction.\n",
91			atypes.first->getName().c_str(),
92			atypes.second->getName().c_str());
93			painCave.severity = OPENMD_ERROR;
94			painCave.isFatal = 1;
95			simError();
96			}
97
98			MorseData* morseData = dynamic_cast<MorseData*>(data);
99			if (morseData == NULL) {
100			sprintf( painCave.errMsg,
101			"Morse::initialize could not convert GenericData to\n"
102			"\tMorseData for %s - %s interaction.\n",
103			atypes.first->getName().c_str(),
104			atypes.second->getName().c_str());
105			painCave.severity = OPENMD_ERROR;
106			painCave.isFatal = 1;
107			simError();
108			}
109
110			MorseParam morseParam = morseData->getData();
111
112			RealType De = morseParam.De;
113			RealType Re = morseParam.Re;
114			RealType beta = morseParam.beta;
115			string interactionType = morseParam.interactionType;
116
117			toUpper(interactionType);
118			map<string, MorseInteractionType>::iterator i;
119			i = stringToEnumMap_.find(interactionType);
120			if (i != stringToEnumMap_.end()) {
121			addExplicitInteraction(atypes.first, atypes.second,
122			De, Re, beta, i->second );
123			} else {
124			sprintf( painCave.errMsg,
125			"Morse::initialize found unknown Morse interaction type\n"
126			"\t(%s) for %s - %s interaction.\n",
127			morseParam.interactionType.c_str(),
128			atypes.first->getName().c_str(),
129			atypes.second->getName().c_str());
130			painCave.severity = OPENMD_ERROR;
131			painCave.isFatal = 1;
132			simError();
133			}
134			}
135			}
136			initialized_ = true;
137			}
138
139			void Morse::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
140			RealType De, RealType Re, RealType beta,
141			MorseInteractionType mit) {
142
143			AtomTypeProperties atp = atype1->getATP();
144			MorseMap.insert( pair<int, AtomType*>(atp.ident, atype1) );
145
146			atp = atype2->getATP();
147			MorseMap.insert( pair<int, AtomType*>(atp.ident, atype2) );
148
149			MorseInteractionData mixer;
150			mixer.De = De;
151			mixer.Re = Re;
152			mixer.beta = beta;
153			mixer.interactionType = mit;
154
155			pair<AtomType, AtomType> key1, key2;
156			key1 = make_pair(atype1, atype2);
157			key2 = make_pair(atype2, atype1);
158
159			MixingMap[key1] = mixer;
160			if (key2 != key1) {
161			MixingMap[key2] = mixer;
162			}
163			}
164
165			void Morse::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
166			RealType r, RealType r2, RealType rcut, RealType sw,
167			RealType vdwMult, RealType &vpair, RealType &pot,
168			Vector3d &f1) {
169
170			if (!initialized_) initialize();
171
172			RealType myPot = 0.0;
173			RealType myPotC = 0.0;
174			RealType myDeriv = 0.0;
175			RealType myDerivC = 0.0;
176
177			pair<AtomType, AtomType> key = make_pair(at1, at2);
178			MorseInteractionData mixer = MixingMap[key];
179
180			RealType De = mixer.De;
181			RealType Re = mixer.Re;
182			RealType beta = mixer.beta;
183			MorseInteractionType interactionType = mixer.interactionType;
184
185			// V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
186
187			RealType expt = -beta*(r - Re);
188			RealType expfnc = exp(expt);
189			RealType expfnc2 = expfnc*expfnc;
190
191			RealType exptC = 0.0;
192			RealType expfncC = 0.0;
193			RealType expfnc2C = 0.0;
194
195			if (Morse::shiftedPot_ \|\| Morse::shiftedFrc_) {
196			exptC = -beta*(rcut - Re);
197			expfncC = exp(exptC);
198			expfnc2C = expfncC*expfncC;
199			}
200
201
202			switch(interactionType) {
203			case shiftedMorse : {
204
205			myPot = De * (expfnc2 - 2.0 * expfnc);
206			myDeriv = 2.0 * De * beta * (expfnc - expfnc2);
207
208			if (Morse::shiftedPot_) {
209			myPotC = De * (expfnc2C - 2.0 * expfncC);
210			myDerivC = 0.0;
211			} else if (Morse::shiftedFrc_) {
212			myPotC = De * (expfnc2C - 2.0 * expfncC);
213			myDerivC = 2.0 * De * beta * (expfnc2C - expfnc2C);
214			myPotC += myDerivC * (r - rcut);
215			} else {
216			myPotC = 0.0;
217			myDerivC = 0.0;
218			}
219
220			break;
221			}
222			case repulsiveMorse : {
223
224			myPot = De * expfnc2;
225			myDeriv = -2.0 * De * beta * expfnc2;
226
227			if (Morse::shiftedPot_) {
228			myPotC = De * expfnc2C;
229			myDerivC = 0.0;
230			} else if (Morse::shiftedFrc_) {
231			myPotC = De * expfnc2C;
232			myDerivC = -2.0 * De * beta * expfnc2C;
233			myPotC += myDerivC * (r - rcut);
234			} else {
235			myPotC = 0.0;
236			myDerivC = 0.0;
237			}
238
239			break;
240			}
241			}
242
243			RealType pot_temp = vdwMult * (myPot - myPotC);
244			vpair += pot_temp;
245
246			RealType dudr = sw * vdwMult * (myDeriv - myDerivC);
247
248			pot += sw * pot_temp;
249			f1 = d * dudr / r;
250
251			return;
252
253			}
254
255			void Morse::do_morse_pair(int atid1, int atid2, RealType *d,
256			RealType rij, RealType r2, RealType *rcut,
257			RealType sw, RealType vdwMult,
258			RealType vpair, RealType pot, RealType *f1) {
259
260			if (!initialized_) initialize();
261
262			AtomType* atype1 = MorseMap[*atid1];
263			AtomType* atype2 = MorseMap[*atid2];
264
265			Vector3d disp(d[0], d[1], d[2]);
266			Vector3d frc(f1[0], f1[1], f1[2]);
267
268			calcForce(atype1, atype2, disp, rij, r2, rcut, sw, vdwMult, vpair,
269			*pot, frc);
270
271			f1[0] = frc.x();
272			f1[1] = frc.y();
273			f1[2] = frc.z();
274
275			return;
276			}
277
278			void Morse::setMorseDefaultCutoff(RealType thisRcut, int sP, int *sF) {
279			shiftedPot_ = (bool)(*sP);
280			shiftedFrc_ = (bool)(*sF);
281			}
282			}
283
284			extern "C" {
285
286			#define fortranSetMorseCutoff FC_FUNC(setmorsedefaultcutoff, SETMORSEDEFAULTCUTOFF)
287			#define fortranDoMorsePair FC_FUNC(do_morse_pair, DO_MORSE_PAIR)
288
289			void fortranSetMorseCutoff(RealType rcut, int shiftedPot, int *shiftedFrc) {
290			return OpenMD::Morse::Instance()->setMorseDefaultCutoff(rcut, shiftedPot,
291			shiftedFrc);
292			}
293			void fortranDoMorsePair(int atid1, int atid2, RealType d, RealType rij,
294			RealType r2, RealType rcut, RealType *sw,
295			RealType vdwMult, RealType vpair, RealType* pot,
296			RealType *f1){
297
298			return OpenMD::Morse::Instance()->do_morse_pair(atid1, atid2, d, rij, r2,
299			rcut, sw, vdwMult, vpair,
300			pot, f1);
301			}
302			}