src/UseTheForce/EAM.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */

#include "UseTheForce/EAM.hpp"

namespace oopse {

//definition of createEAM
ForceField* createEAM() {
    return new EAM();
}

//register createDUFF to ForceFieldFactory
ForceFieldFactory::getInstance()->registerForceField("EAM", createDUFF);

void EAM::parse(const std::string& filename) {
    ifstrstream* ffStream;
    ffStream = openForceFiledFile(filename);
    const int bufferSize = 65535;
    char buffer[bufferSize];
    std::string line;

    int status;
    int ident = 1; //fortran's index begins from 1
    int lineNo = 0;

    while(ffStream.getline(buffer, bufferSize)){
        ++lineNo;
        line = trimSpaces(buffer);
        
        if ( line.empty() || (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
            continue;
        } else {
            AtomType* atomType;
            std::string atomTypeName;
            double mass;        
            StringTokenizer tokenizer(line);
            
            if (tokenizer.countToken() >= 3) {
                atomTypeName = tokenizer.nextToken();
                mass = tokenizer.nextTokenAsDouble();
                std::string potentialParamFile = tokenizer.nextToken();
                
                atomType = new AtomType();
                atomType->setName(atomTypeName);
                atomType->setMass(mass);                

                atomType->setIdent(ident);
                atomType->setEAM();                            
                atomType->complete();

                parseEAMParamFile(potentialParamFile);
                                
                //add atom type to AtomTypeContainer
                addAtomType(atomTypeName, atomType);
                ++ident;

                
            } else {
                sprintf( painCave.errMsg,
                       "Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
                       "in line %d : %s\n",
                       filename.c_str(), lineNo, line);
                painCave.severity = OOPSE_ERROR;
                painCave.isFatal = 1;
                simError();                
            }
            
        }

    }
    
    delete ffStream;

}

void EAM::parseEAMParamFile(const std::string& potentialParamFile) {

    ifstrstream* ppfStream;
    ppfStream = openForceFiledFile(potentialParamFile);
    const int bufferSize = 65535;
    char buffer[bufferSize];
    std::string line;

    //skip first line
    ppfStream->getline(buffer, bufferSize);


    //The Second line contains atomic number, atomic mass, a lattice constant and lattic type
    int ident;
    double mass;
    double latticeConstant; 
    std::string lattic;
    if (ppfStream->getline(buffer, bufferSize)) {       
        StringTokenizer tokenizer1(buffer);
        ident = tokenizer1.nextTokenAsInt();
        
        if (tokenizer1.countTokens() >= 4) {
            ident = tokenizer1.nextTokenAsInt();
            mass = tokenizer1.nextTokenAsDouble();
            latticeConstant = tokenizer1.nextTokenAsDouble();
            lattice = tokenizer1.nextToken();
        }else {
            std::cerr << "Not enought tokens" << std::endl;
        }
    } else {

    }
    
    // The third line is nrho, drho, nr, dr and rcut
    int nrho;
    double drho;
    int nr;
    double dr;
    double rcut;
    
    if (ppfStream->getline(buffer, bufferSize)) {
        StringTokenizer tokenizer2(buffer);

        if (tokenizer2.countTokens() >= 5){
            nrho = tokenizer2.nextTokenAsInt();
            drho = tokenizer2.nextTokenAsDouble();
            nr = tokenizer2.nextTokenAsInt();
            dr = tokenizer2.nextTokenAsDouble();
            rcut = tokenizer2.nextTokenAsDouble();
        }else {
            std::cerr << "Not enought tokens" << std::endl;
        }
    } else {

    }
    
    std::vector<double> rvals; 
    std::vector<double> rhovals;
    std::vector<double> Frhovals;    

    parseEAMArray(ppfStream, rvals, nr);
    parseEAMArray(ppfStream, rhovals, nr);
    parseEAMArray(ppfStream, Frhovals, nrho);

    int status;
    newEAMtype(&latticeConstant, &nrho, &drho,  &nr, &dr, &rcut,
                        &rvals[0], &rhovals[0], &Frhovals[0], &ident, &status );

    if (status == 0) {

    }
       
}

void EAM::parseEAMArray(ifstrstream* ppfStream, std::vector<double>& array, int num) {
    
    const int dataPerLine = 5;
    if (num % dataPerLine != 0) {

    }
    
    const int bufferSize = 65535;
    char buffer[bufferSize];
    std::string line;
    int readLines = num/dataPerLine;
    int lineCount = 0;
    while(ffStream.getline(buffer, bufferSize) && lineCount < num){

        StringTokenizer tokenizer(buffer);
        if (tokenizer.countTokens() >= 5) {
            array.push_back(tokenizer.nextTokenAsDouble));
        } else {

        }
        ++lineCount;
    }

    if (lineCount < num) {
        
    }
    
}

} //end namespace oopse
Revision:	1760
Committed:	Fri Nov 19 20:23:26 2004 UTC (19 years, 9 months ago) by tim
File size:	6049 byte(s)
Log Message:	add EAM
#	User	Rev	Content
1	tim	1760	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.org
5			*
6			* This program is free software; you can redistribute it and/or
7			* modify it under the terms of the GNU Lesser General Public License
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
22			* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23			*
24			*/
25
26			#include "UseTheForce/EAM.hpp"
27
28			namespace oopse {
29
30			//definition of createEAM
31			ForceField* createEAM() {
32			return new EAM();
33			}
34
35			//register createDUFF to ForceFieldFactory
36			ForceFieldFactory::getInstance()->registerForceField("EAM", createDUFF);
37
38			void EAM::parse(const std::string& filename) {
39			ifstrstream* ffStream;
40			ffStream = openForceFiledFile(filename);
41			const int bufferSize = 65535;
42			char buffer[bufferSize];
43			std::string line;
44
45			int status;
46			int ident = 1; //fortran's index begins from 1
47			int lineNo = 0;
48
49			while(ffStream.getline(buffer, bufferSize)){
50			++lineNo;
51			line = trimSpaces(buffer);
52
53			if ( line.empty() \|\| (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
54			continue;
55			} else {
56			AtomType* atomType;
57			std::string atomTypeName;
58			double mass;
59			StringTokenizer tokenizer(line);
60
61			if (tokenizer.countToken() >= 3) {
62			atomTypeName = tokenizer.nextToken();
63			mass = tokenizer.nextTokenAsDouble();
64			std::string potentialParamFile = tokenizer.nextToken();
65
66			atomType = new AtomType();
67			atomType->setName(atomTypeName);
68			atomType->setMass(mass);
69
70			atomType->setIdent(ident);
71			atomType->setEAM();
72			atomType->complete();
73
74			parseEAMParamFile(potentialParamFile);
75
76			//add atom type to AtomTypeContainer
77			addAtomType(atomTypeName, atomType);
78			++ident;
79
80
81			} else {
82			sprintf( painCave.errMsg,
83			"Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
84			"in line %d : %s\n",
85			filename.c_str(), lineNo, line);
86			painCave.severity = OOPSE_ERROR;
87			painCave.isFatal = 1;
88			simError();
89			}
90
91			}
92
93			}
94
95			delete ffStream;
96
97			}
98
99			void EAM::parseEAMParamFile(const std::string& potentialParamFile) {
100
101			ifstrstream* ppfStream;
102			ppfStream = openForceFiledFile(potentialParamFile);
103			const int bufferSize = 65535;
104			char buffer[bufferSize];
105			std::string line;
106
107			//skip first line
108			ppfStream->getline(buffer, bufferSize);
109
110
111			//The Second line contains atomic number, atomic mass, a lattice constant and lattic type
112			int ident;
113			double mass;
114			double latticeConstant;
115			std::string lattic;
116			if (ppfStream->getline(buffer, bufferSize)) {
117			StringTokenizer tokenizer1(buffer);
118			ident = tokenizer1.nextTokenAsInt();
119
120			if (tokenizer1.countTokens() >= 4) {
121			ident = tokenizer1.nextTokenAsInt();
122			mass = tokenizer1.nextTokenAsDouble();
123			latticeConstant = tokenizer1.nextTokenAsDouble();
124			lattice = tokenizer1.nextToken();
125			}else {
126			std::cerr << "Not enought tokens" << std::endl;
127			}
128			} else {
129
130			}
131
132			// The third line is nrho, drho, nr, dr and rcut
133			int nrho;
134			double drho;
135			int nr;
136			double dr;
137			double rcut;
138
139			if (ppfStream->getline(buffer, bufferSize)) {
140			StringTokenizer tokenizer2(buffer);
141
142			if (tokenizer2.countTokens() >= 5){
143			nrho = tokenizer2.nextTokenAsInt();
144			drho = tokenizer2.nextTokenAsDouble();
145			nr = tokenizer2.nextTokenAsInt();
146			dr = tokenizer2.nextTokenAsDouble();
147			rcut = tokenizer2.nextTokenAsDouble();
148			}else {
149			std::cerr << "Not enought tokens" << std::endl;
150			}
151			} else {
152
153			}
154
155			std::vector<double> rvals;
156			std::vector<double> rhovals;
157			std::vector<double> Frhovals;
158
159			parseEAMArray(ppfStream, rvals, nr);
160			parseEAMArray(ppfStream, rhovals, nr);
161			parseEAMArray(ppfStream, Frhovals, nrho);
162
163			int status;
164			newEAMtype(&latticeConstant, &nrho, &drho, &nr, &dr, &rcut,
165			&rvals[0], &rhovals[0], &Frhovals[0], &ident, &status );
166
167			if (status == 0) {
168
169			}
170
171			}
172
173			void EAM::parseEAMArray(ifstrstream* ppfStream, std::vector<double>& array, int num) {
174
175			const int dataPerLine = 5;
176			if (num % dataPerLine != 0) {
177
178			}
179
180			const int bufferSize = 65535;
181			char buffer[bufferSize];
182			std::string line;
183			int readLines = num/dataPerLine;
184			int lineCount = 0;
185			while(ffStream.getline(buffer, bufferSize) && lineCount < num){
186
187			StringTokenizer tokenizer(buffer);
188			if (tokenizer.countTokens() >= 5) {
189			array.push_back(tokenizer.nextTokenAsDouble));
190			} else {
191
192			}
193			++lineCount;
194			}
195
196			if (lineCount < num) {
197
198			}
199
200			}
201
202			} //end namespace oopse