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.countTokens() >= 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:	1762
Committed:	Fri Nov 19 21:38:22 2004 UTC (19 years, 7 months ago) by tim
File size:	6050 byte(s)
Log Message:	refactory integrator
#	Content
1	/*
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.countTokens() >= 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