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, 9 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.
#	Content
1	/*
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	* 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
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	#include "io/BaseAtomTypesSectionParser.hpp"
51	#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	namespace OpenMD {
69
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	spMan_.push_back(new BaseAtomTypesSectionParser());
90	spMan_.push_back(new AtomTypesSectionParser());
91	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	}
106
107	void MnM_FF::parse(const std::string& filename) {
108	ifstrstream* ffStream;
109
110	ffStream = openForceFieldFile(filename);
111
112	spMan_.parse(ffStream, this);
113
114	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	// 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	for (int j = ayb.size()-1; j > 0; j--) {
129	ayb[j-1]->useBase(ayb[j]);
130	}
131	}
132	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	painCave.severity = OPENMD_ERROR;
158	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
168	/* 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
181	/* get the name and ident of the metallic atom type */
182	std::string at1s = at->getName();
183	int atid1 = at->getIdent();
184
185	/* get the name and ident of the nonmetallic atom type */
186	std::string at2s = at2->getName();
187	int atid2 = at2->getIdent();
188
189	/* look for a match in the non-bonded interactions parsed
190	from the force field file */
191	nbit = getNonBondedInteractionType(at1s, at2s);
192
193	/* if we found a match (even a partial match), punt to the
194	interaction to poke our info down to fortran. */
195	if (nbit != NULL) nbit->tellFortran(atid1, atid2);
196	}
197	}
198	}
199	}
200
201	delete ffStream;
202
203	}
204
205
206	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	painCave.severity = OPENMD_ERROR;
222	painCave.isFatal = 1;
223	simError();
224	}
225	}
226	else {
227	sprintf(painCave.errMsg, "Can not find EAM Parameters\n");
228	painCave.severity = OPENMD_ERROR;
229	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	} //end namespace OpenMD
247