src/UseTheForce/MnM_FF.cpp

/*
 * Copyright (c) 2007 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 "UseTheForce/MnM_FF.hpp"
#include "UseTheForce/DarkSide/lj_interface.h"
#include "UseTheForce/DarkSide/sticky_interface.h"
#include "UseTheForce/DarkSide/eam_interface.h"
#include "UseTheForce/DarkSide/suttonchen_interface.h"
#include "UseTheForce/ForceFieldFactory.hpp"
#include "io/DirectionalAtomTypesSectionParser.hpp"
#include "io/BaseAtomTypesSectionParser.hpp"
#include "io/AtomTypesSectionParser.hpp"
#include "io/LennardJonesAtomTypesSectionParser.hpp"
#include "io/ChargeAtomTypesSectionParser.hpp"
#include "io/MultipoleAtomTypesSectionParser.hpp"
#include "io/StickyAtomTypesSectionParser.hpp"
#include "io/StickyPowerAtomTypesSectionParser.hpp"
#include "io/GayBerneAtomTypesSectionParser.hpp"
#include "io/BondTypesSectionParser.hpp"
#include "io/BendTypesSectionParser.hpp"
#include "io/TorsionTypesSectionParser.hpp"
#include "io/MetalNonMetalInteractionsSectionParser.hpp"
#include "io/EAMAtomTypesSectionParser.hpp"
#include "io/SCAtomTypesSectionParser.hpp"
#include "io/OptionSectionParser.hpp"
#include "UseTheForce/ForceFieldCreator.hpp"


namespace OpenMD {
   
  MnM_FF::MnM_FF() {    
    
    //set default force field filename
    setForceFieldFileName("MnM.frc");
    
    //The order of adding section parsers is important.
    //OptionSectionParser must come first to set options for other parsers
    //DirectionalAtomTypesSectionParser should be added before
    //AtomTypesSectionParser, and these two section parsers will actually
    //create "real" AtomTypes (AtomTypesSectionParser will create AtomType and
    //DirectionalAtomTypesSectionParser will create DirectionalAtomType, which
    //is a subclass of AtomType and should come first). Other AtomTypes Section
    //Parser will not create the "real" AtomType, they only add and set some
    //attribute of the AtomType. Thus their order are not important.
    //AtomTypesSectionParser should be added before other atom type section
    //parsers. Make sure they are added after DirectionalAtomTypesSectionParser
    //and AtomTypesSectionParser. The order of BondTypesSectionParser,
    //BendTypesSectionParser and TorsionTypesSectionParser are not important.
    spMan_.push_back(new OptionSectionParser(forceFieldOptions_));
    spMan_.push_back(new BaseAtomTypesSectionParser());
    spMan_.push_back(new AtomTypesSectionParser());
    spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new LennardJonesAtomTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new MetalNonMetalInteractionsSectionParser(forceFieldOptions_));
    spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_));
    spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_));
    
  }
  
  void MnM_FF::parse(const std::string& filename) {
    ifstrstream* ffStream;
        
    ffStream = openForceFieldFile(filename);
    
    spMan_.parse(*ffStream, *this);
    
    ForceField::AtomTypeContainer::MapTypeIterator i;
    AtomType* at;
    ForceField::AtomTypeContainer::MapTypeIterator j;
    AtomType* at2;
    ForceField::NonBondedInteractionTypeContainer::MapTypeIterator k;
    NonBondedInteractionType* nbit;
    
    for (at = atomTypeCont_.beginType(i); at != NULL;
         at = atomTypeCont_.nextType(i)) {
      // useBase sets the responsibilities, and these have to be done 
      // after the atomTypes and Base types have all been scanned:

      std::vector<AtomType*> ayb = at->allYourBase();      
      if (ayb.size() > 1) {
        for (int j = ayb.size()-1; j > 0; j--) {          
          ayb[j-1]->useBase(ayb[j]);
        }
      }
      at->makeFortranAtomType();
    }
    
    for (at = atomTypeCont_.beginType(i); at != NULL;
         at = atomTypeCont_.nextType(i)) {
      at->complete();
    }
    
    hasSCtypes_ = false;
    for (at = atomTypeCont_.beginType(i); at != NULL;
         at = atomTypeCont_.nextType(i)) {
      if (at->isSC())
        hasSCtypes_ = true;
    }
    
    hasEAMtypes_ = false;
    for (at = atomTypeCont_.beginType(i); at != NULL;
         at = atomTypeCont_.nextType(i)) {
      if (at->isEAM())
        hasEAMtypes_ = true;
    }
    
    if (hasEAMtypes_ && hasSCtypes_) {
      sprintf(painCave.errMsg,
              "MnM_FF forcefield cannot use both EAM and Sutton-Chen at the same time\n");
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();
    }
    
    /* to handle metal-nonmetal interactions, first we loop over
       all atom types: */ 
    
    for (at = atomTypeCont_.beginType(i); at != NULL;
         at = atomTypeCont_.nextType(i)) {
      
      /* if we find a metallic atom, we need to compare against
         all other atom types */ 
      
      if (at->isEAM() || at->isSC()) {
        
        /* loop over all other atom types */
        for (at2 = atomTypeCont_.beginType(j); at2 != NULL;
             at2 = atomTypeCont_.nextType(j)) {
          
          /* if the other partner is not a metallic type, we need to
             look for explicit non-bonded interactions */
          if (!at2->isEAM() && !at2->isSC()) {
          
            /* get the name and ident of the metallic atom type */
            std::string at1s = at->getName();
            int atid1 = at->getIdent();
          
            /* get the name and ident of the nonmetallic atom type */
            std::string at2s = at2->getName();
            int atid2 = at2->getIdent();
          
            /* look for a match in the non-bonded interactions parsed
               from the force field file */ 
            nbit = getNonBondedInteractionType(at1s, at2s);
          
            /* if we found a match (even a partial match), punt to the
               interaction to poke our info down to fortran. */
            if (nbit != NULL) nbit->tellFortran(atid1, atid2);
          }
        }                       
      }
    }
    
    delete ffStream;
    
  }
  
  
  RealType MnM_FF::getRcutFromAtomType(AtomType* at) {
    RealType rcut = 0.0;
    if (at->isEAM()) {
      GenericData* data = at->getPropertyByName("EAM");
      if (data != NULL) {
        EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
        
        if (eamData != NULL) {
          
          EAMParam& eamParam = eamData->getData();
          rcut =  eamParam.rcut;
        }
        else {
          sprintf(painCave.errMsg,
                  "Can not cast GenericData to EAMParam\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError();
        }
      }
      else {
        sprintf(painCave.errMsg, "Can not find EAM Parameters\n");
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();
      }
    }
    else {
      rcut = ForceField::getRcutFromAtomType(at);
    }
    
    return rcut;
  }
     
  MnM_FF::~MnM_FF() {
    destroyLJTypes();
    destroyStickyTypes();
    if (hasEAMtypes_) destroyEAMTypes();
    if (hasSCtypes_)  destroySCTypes();
  }
} //end namespace OpenMD

Revision:	1596
Committed:	Mon Jul 25 17:30:53 2011 UTC (13 years, 10 months ago) by gezelter
File size:	9721 byte(s)
Log Message:	Updated the BlockSnapshotManager to use a specified memory footprint in constructor and not to rely on physmem and residentMem to figure out free memory. DynamicProps is the only program that uses the BlockSnapshotManager, so substantial changes were needed there as well.
#	User	Rev	Content
1	chuckv	1152	/*
2			* Copyright (c) 2007 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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	chuckv	1152	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	chuckv	1152	* 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	gezelter	1390	*
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	chuckv	1152	*/
41
42
43			#include "UseTheForce/MnM_FF.hpp"
44			#include "UseTheForce/DarkSide/lj_interface.h"
45			#include "UseTheForce/DarkSide/sticky_interface.h"
46			#include "UseTheForce/DarkSide/eam_interface.h"
47			#include "UseTheForce/DarkSide/suttonchen_interface.h"
48			#include "UseTheForce/ForceFieldFactory.hpp"
49			#include "io/DirectionalAtomTypesSectionParser.hpp"
50	gezelter	1280	#include "io/BaseAtomTypesSectionParser.hpp"
51	chuckv	1152	#include "io/AtomTypesSectionParser.hpp"
52			#include "io/LennardJonesAtomTypesSectionParser.hpp"
53			#include "io/ChargeAtomTypesSectionParser.hpp"
54			#include "io/MultipoleAtomTypesSectionParser.hpp"
55			#include "io/StickyAtomTypesSectionParser.hpp"
56			#include "io/StickyPowerAtomTypesSectionParser.hpp"
57			#include "io/GayBerneAtomTypesSectionParser.hpp"
58			#include "io/BondTypesSectionParser.hpp"
59			#include "io/BendTypesSectionParser.hpp"
60			#include "io/TorsionTypesSectionParser.hpp"
61			#include "io/MetalNonMetalInteractionsSectionParser.hpp"
62			#include "io/EAMAtomTypesSectionParser.hpp"
63			#include "io/SCAtomTypesSectionParser.hpp"
64			#include "io/OptionSectionParser.hpp"
65			#include "UseTheForce/ForceFieldCreator.hpp"
66
67
68	gezelter	1390	namespace OpenMD {
69	gezelter	1238
70			MnM_FF::MnM_FF() {
71
72			//set default force field filename
73			setForceFieldFileName("MnM.frc");
74
75			//The order of adding section parsers is important.
76			//OptionSectionParser must come first to set options for other parsers
77			//DirectionalAtomTypesSectionParser should be added before
78			//AtomTypesSectionParser, and these two section parsers will actually
79			//create "real" AtomTypes (AtomTypesSectionParser will create AtomType and
80			//DirectionalAtomTypesSectionParser will create DirectionalAtomType, which
81			//is a subclass of AtomType and should come first). Other AtomTypes Section
82			//Parser will not create the "real" AtomType, they only add and set some
83			//attribute of the AtomType. Thus their order are not important.
84			//AtomTypesSectionParser should be added before other atom type section
85			//parsers. Make sure they are added after DirectionalAtomTypesSectionParser
86			//and AtomTypesSectionParser. The order of BondTypesSectionParser,
87			//BendTypesSectionParser and TorsionTypesSectionParser are not important.
88			spMan_.push_back(new OptionSectionParser(forceFieldOptions_));
89	gezelter	1280	spMan_.push_back(new BaseAtomTypesSectionParser());
90	gezelter	1282	spMan_.push_back(new AtomTypesSectionParser());
91	gezelter	1238	spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_));
92			spMan_.push_back(new LennardJonesAtomTypesSectionParser(forceFieldOptions_));
93			spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_));
94			spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_));
95			spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_));
96			spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_));
97			spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_));
98			spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_));
99			spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_));
100			spMan_.push_back(new MetalNonMetalInteractionsSectionParser(forceFieldOptions_));
101			spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_));
102			spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_));
103			spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_));
104
105	chuckv	1162	}
106
107	gezelter	1238	void MnM_FF::parse(const std::string& filename) {
108			ifstrstream* ffStream;
109
110			ffStream = openForceFieldFile(filename);
111	chuckv	1152
112	gezelter	1238	spMan_.parse(ffStream, this);
113	chuckv	1162
114	gezelter	1238	ForceField::AtomTypeContainer::MapTypeIterator i;
115			AtomType* at;
116			ForceField::AtomTypeContainer::MapTypeIterator j;
117			AtomType* at2;
118			ForceField::NonBondedInteractionTypeContainer::MapTypeIterator k;
119			NonBondedInteractionType* nbit;
120
121			for (at = atomTypeCont_.beginType(i); at != NULL;
122			at = atomTypeCont_.nextType(i)) {
123	gezelter	1280	// useBase sets the responsibilities, and these have to be done
124			// after the atomTypes and Base types have all been scanned:
125
126			std::vector<AtomType*> ayb = at->allYourBase();
127			if (ayb.size() > 1) {
128	gezelter	1555	for (int j = ayb.size()-1; j > 0; j--) {
129	gezelter	1280	ayb[j-1]->useBase(ayb[j]);
130			}
131			}
132	gezelter	1238	at->makeFortranAtomType();
133			}
134
135			for (at = atomTypeCont_.beginType(i); at != NULL;
136			at = atomTypeCont_.nextType(i)) {
137			at->complete();
138			}
139
140			hasSCtypes_ = false;
141			for (at = atomTypeCont_.beginType(i); at != NULL;
142			at = atomTypeCont_.nextType(i)) {
143			if (at->isSC())
144			hasSCtypes_ = true;
145			}
146
147			hasEAMtypes_ = false;
148			for (at = atomTypeCont_.beginType(i); at != NULL;
149			at = atomTypeCont_.nextType(i)) {
150			if (at->isEAM())
151			hasEAMtypes_ = true;
152			}
153
154			if (hasEAMtypes_ && hasSCtypes_) {
155			sprintf(painCave.errMsg,
156			"MnM_FF forcefield cannot use both EAM and Sutton-Chen at the same time\n");
157	gezelter	1390	painCave.severity = OPENMD_ERROR;
158	gezelter	1238	painCave.isFatal = 1;
159			simError();
160			}
161
162			/* to handle metal-nonmetal interactions, first we loop over
163			all atom types: */
164
165			for (at = atomTypeCont_.beginType(i); at != NULL;
166			at = atomTypeCont_.nextType(i)) {
167	chuckv	1162
168	gezelter	1238	/* if we find a metallic atom, we need to compare against
169			all other atom types */
170
171			if (at->isEAM() \|\| at->isSC()) {
172
173			/* loop over all other atom types */
174			for (at2 = atomTypeCont_.beginType(j); at2 != NULL;
175			at2 = atomTypeCont_.nextType(j)) {
176
177			/* if the other partner is not a metallic type, we need to
178			look for explicit non-bonded interactions */
179			if (!at2->isEAM() && !at2->isSC()) {
180	chuckv	1162
181	gezelter	1238	/* get the name and ident of the metallic atom type */
182			std::string at1s = at->getName();
183			int atid1 = at->getIdent();
184	chuckv	1162
185	gezelter	1238	/* get the name and ident of the nonmetallic atom type */
186			std::string at2s = at2->getName();
187			int atid2 = at2->getIdent();
188	chuckv	1162
189	gezelter	1238	/* look for a match in the non-bonded interactions parsed
190			from the force field file */
191			nbit = getNonBondedInteractionType(at1s, at2s);
192	chuckv	1252
193	gezelter	1238	/* if we found a match (even a partial match), punt to the
194			interaction to poke our info down to fortran. */
195	gezelter	1596	if (nbit != NULL) nbit->tellFortran(atid1, atid2);
196			}
197			}
198	chuckv	1162	}
199	chuckv	1152	}
200	gezelter	1596
201	gezelter	1238	delete ffStream;
202	gezelter	1596
203	chuckv	1152	}
204	chuckv	1162
205
206	gezelter	1238	RealType MnM_FF::getRcutFromAtomType(AtomType* at) {
207			RealType rcut = 0.0;
208			if (at->isEAM()) {
209			GenericData* data = at->getPropertyByName("EAM");
210			if (data != NULL) {
211			EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
212
213			if (eamData != NULL) {
214
215			EAMParam& eamParam = eamData->getData();
216			rcut = eamParam.rcut;
217			}
218			else {
219			sprintf(painCave.errMsg,
220			"Can not cast GenericData to EAMParam\n");
221	gezelter	1390	painCave.severity = OPENMD_ERROR;
222	gezelter	1238	painCave.isFatal = 1;
223			simError();
224			}
225			}
226			else {
227			sprintf(painCave.errMsg, "Can not find EAM Parameters\n");
228	gezelter	1390	painCave.severity = OPENMD_ERROR;
229	gezelter	1238	painCave.isFatal = 1;
230			simError();
231			}
232			}
233			else {
234			rcut = ForceField::getRcutFromAtomType(at);
235			}
236
237			return rcut;
238			}
239
240			MnM_FF::~MnM_FF() {
241			destroyLJTypes();
242			destroyStickyTypes();
243			if (hasEAMtypes_) destroyEAMTypes();
244			if (hasSCtypes_) destroySCTypes();
245			}
246	gezelter	1390	} //end namespace OpenMD
247	chuckv	1162