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 (15 years, 2 months ago) by gezelter
Original Path:	*branches/development/src/nonbonded/Morse.cpp*
File size:	10196 byte(s)
Log Message:	Starting migration of Morse to C++
#	Content
1	/*
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	}