src/UseTheForce/DUFF.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/DUFF.hpp"
#include "UseTheForce/DarkSide/lj_interface.h"

namespace oopse {

//definition of createDUFF
ForceField* createDUFF() {
    return new DUFF();
}

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


ParseState DUFF::getSection(const std::string& section) {
    ParseState result;
    
    switch(section) {
        case "AtomTypes" :
            result = DUFF::AtomTypeSection;
            break;

        case "BondTypes" :
            result = DUFF::BondTypeSection;
            break;

        case "BendTypes" :
            result = DUFF::BendTypeSection;
            break;

        case "TorsionTypes" :
            result = DUFF::TorsionTypeSection;
            break;
        default:
            result = DUFF::UnknownSection;
    }

    return result;
}

void DUFF::parse(const std::string& filename) {
    ifstrstream* ffStream;
    ffStream = openForceFiledFile(filename);
    const int bufferSize = 65535;
    std::string line;
    char buffer[bufferSize];
    int lineNo = 0;
    ParseState currentSection = DUFF::UnknownSection;
    
    while(ffStream.getline(buffer, bufferSize)){
        ++lineNo;

        line = trimSpaces(buffer);
        //a line begins with "//" is comment
        if ( line.empty() || (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
            continue;
        } else {

            switch(currentSection) {
                case DUFF::AtomTypeSection :
                    parseAtomType(line, lineNo);
                    break;

                case DUFF::BondTypeSection :
                    parseBondType(line, lineNo);
                    break;

                case DUFF::BendTypeSection :
                    parseBendType(line, lineNo);
                    break;

                case DUFF::TorsionTypeSection :
                    parseTorsionType(line, lineNo);
                    break;
                    
                case DUFF::UnknownSection:
                    StringTokenizer tokenizer(line);
                    
                    std::string keyword = tokenizer.nextToken();
                    std::string section = tokenizer.nextToken();

                    ParseState newSection = getSection(section);
                    if (keyword != "begin" || keyword != "end") {
                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
                    } else if (keyword == "begin") {
                        if (newSection == DUFF::UnknownSection) {
                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
                        } else {
                            //enter a new section
                            currentSection = newSection;
                        }
                        
                    } else if (keyword == "end"){
                        if (currentSection == newSection) ) {
                            //leave a section
                            currentSection = DUFF::UnknownSection;
                        } else {                            
                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
                        }

                    }
                    break;
                default :
                
            }
            
        }
    }
    
    delete ffStream;
}


void DUFF::parseAtomType(const std::string& line, int lineNo){ 
    StringTokenizer tokenizer(line);
    
    if (tokenizer.countTokens() >= 4) {
        atomTypeName = tokenizer.nextToken();
        mass = tokenizer.nextTokenAsFloat();
        epsilon = tokenizer.nextTokenAsFloat();
        sigma = tokenizer.nextTokenAsFloat();

        atomType = new AtomType();
        atomType->setName(atomTypeName);
        atomType->setMass(mass);
        
        //by default, all of the properties in AtomTypeProperties is set to 0
        //In  Lennard-Jones Force Field, we only need to set Lennard-Jones to true
        atomType->setLennardJones();

        atomType->setIdent(ident);
        atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
        atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
        atomType->complete();
        //notify a new LJtype atom type is created
        newLJtype(&ident, &sigma, &epsilon, &status);

        //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();                
    }    
}

void DUFF::parseBondType(const std::string& line, int lineNo){ 

    StringTokenizer tokenizer(line);
    std::string at1;
    std::string at2;
    std::string bt;
    BondType* bondType = NULL;
    double b0;
    
    int nTokens = tokenizer.countTokens();

    if (nTokens < 4) {

        return;
    }
    
    at1 = tokenizer.nextToken();
    at2 = tokenizer.nextToken();
    bt = tokenizer.nextToken();
    b0 = tokenizer.nextTokenAsDouble();
    nTokens -= 4;

    //switch is a maintain nightmare
    switch(bt) {
        case "FixedBondType" :
            bondType = new FixedBondType();
            break;
            
        case "HarmonicBondType" :
            if (nTokens < 1) {

            } else {

                double kb = tokenizer.nextTokenAsDouble();
                bondType = new HarmonicBondType(b0, kb);
            }

            break;

        case "CubicBondType" :
            if (nTokens < 4) {

            } else {

                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bondType = new CubicBondType(b0, k3, k2, k1, k0);
            }
            break;
            
        case "QuadraticBondType" :
            if (nTokens < 5) {

            } else {

                b0 = tokenizer.nextTokenAsDouble();
                double k4 = tokenizer.nextTokenAsDouble();
                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bondType = new QuadraticBondType(b0, k4, k3, k2, k1, k0);
            }
            break;

        case "PolynomialBondType " :
            if (nTokens < 2 || nTokens % 2 != 0) {

            } else {
                int nPairs = nTokens / 2;
                int power;
                double coefficient;
                PolynomialBondType pbt = new PolynomialBondType();
                
                for (int i = 0; i < nPairs; ++i) {
                    power = tokenizer.nextTokenAsInt();
                    coefficient = tokenizer.nextTokenAsDouble();

                    if (power < 0) {

                    } else {
                        pbt->setCoefficient(power, coefficient);
                    }
                }
            }
            
            break;

        default:
            
    }

    if (bondType != NULL) {
        addBondType(at1, at2, bondType);
    }
}

void DUFF::parseBendType(const std::string& line, int lineNo){ 
    StringTokenizer tokenizer(line);
    std::string at1;
    std::string at2;
    std::string at3;
    std::string bt;
    double theta0;
    BendType* bendType = NULL;

    int nTokens = tokenizer.countTokens();

    if (nTokens < 5) {

        return;
    }
    
    at1 = tokenizer.nextToken();
    at2 = tokenizer.nextToken();
    at3 = tokenizer.nextToken();
    bt = tokenizer.nextToken();
    theta0 = tokenizer.nextTokenAsDouble();
    nTokens -= 5;

    //switch is a maintain nightmare
    switch(bt) {
            
        case "HarmonicBendType" :
            
            if (nTokens < 1) {

            } else {

                double ktheta = tokenizer.nextTokenAsDouble();
                bendType = new HarmonicBendType(theta0, ktheta);
            }
            break;
        case "GhostBendType" :
            if (nTokens < 1) {

            } else {
                double ktheta = tokenizer.nextTokenAsDouble();
                bendType = new HarmonicBendType(theta0, ktheta);                
            }
            break;            

        case "UreyBradleyBendType" :
            if (nTokens < 3) {

            } else {
                double ktheta = tokenizer.nextTokenAsDouble();
                double s0 =  tokenizer.nextTokenAsDouble();
                double kub = tokenizer.nextTokenAsDouble();
                bendType = new UreyBradleyBendType(theta0, ktheta, s0, kub);                
            }
            break; 
            
        case "CubicBendType" :
            if (nTokens < 4) {

            } else {

                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bendType = new CubicBendType(theta0, k3, k2, k1, k0);
            }
            break;
            
        case "QuadraticBendType" :
            if (nTokens < 5) {

            } else {

                theta0 = tokenizer.nextTokenAsDouble();
                double k4 = tokenizer.nextTokenAsDouble();
                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bendType = new QuadraticBendType(theta0, k4, k3, k2, k1, k0);
            }
            break;

        case "PolynomialBendType " :
            if (nTokens < 2 || nTokens % 2 != 0) {

            } else {
                int nPairs = nTokens / 2;
                int power;
                double coefficient;
                PolynomialBendType* pbt = new PolynomialBendType();
                
                for (int i = 0; i < nPairs; ++i) {
                    power = tokenizer.nextTokenAsInt();
                    coefficient = tokenizer.nextTokenAsDouble();

                    if (power < 0) {

                    } else {
                        pbt->setCoefficient(power, coefficient);
                    }
                }
            }
            
            break;

        default:
            
    }

    if (bendType != NULL) {
        addBendType(at1, at2, at3, bendType);
    }
    
}

void DUFF::parseTorsionType(const std::string& line, int lineNo){ 
    StringTokenizer tokenizer(line);
    std::string at1;
    std::string at2;
    std::string at3;
    std::string at4;
    std::string tt;
    TorsionType* bendType = NULL;

    int nTokens = tokenizer.countTokens();

    if (nTokens < 5) {

        return;
    }
    
    at1 = tokenizer.nextToken();
    at2 = tokenizer.nextToken();
    at3 = tokenizer.nextToken();
    at4 = tokenizer.nextToken();
    tt = tokenizer.nextToken();

    nTokens -= 5;

    switch(tt) {
            
        case "CubicTorsionType" :
            if (nTokens < 4) {

            } else {

                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bendType = new CubicTorsionType(k3, k2, k1, k0);
            }
            break;
            
        case "QuadraticTorsionType" :
            if (nTokens < 5) {

            } else {

                theta0 = tokenizer.nextTokenAsDouble();
                double k4 = tokenizer.nextTokenAsDouble();
                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bendType = new QuadraticTorsionType( k4, k3, k2, k1, k0);
            }
            break;

        case "PolynomialTorsionType " :
            if (nTokens < 2 || nTokens % 2 != 0) {

            } else {
                int nPairs = nTokens / 2;
                int power;
                double coefficient;
                PolynomialTorsionType* pbt = new PolynomialTorsionType();
                
                for (int i = 0; i < nPairs; ++i) {
                    power = tokenizer.nextTokenAsInt();
                    coefficient = tokenizer.nextTokenAsDouble();

                    if (power < 0) {

                    } else {
                        pbt->setCoefficient(power, coefficient);
                    }
                }
            }
            
            break;
        case "CharmmTorsionType" :
            
            if (nTokens < 3 || nTokens % 3 != 0) {

            } else {
                int nSets = nTokens / 3;
  
                double kchi;
                int n;
                double delta;                
                CharmmTorsionType* ctt = new CharmmTorsionType();
                
                for (int i = 0; i < nSets; ++i) {
                    kchi = tokenizer.nextTokenAsDouble();
                    n = tokenizer.nextTokenAsInt();
                    delta = tokenizer.nextTokenAsDouble();
    
                    if (n < 0 || n > 6) {

                    } else {
                        ctt->setCharmmTorsionParameter(kchi, n, delta);
                    }
                }
            }
        default:
            
    }

    if (bendType != NULL) {
        addTorsionType(at1, at2, at3, bendType);
    }
}

} //end namespace oopse

Revision:	1757
Committed:	Fri Nov 19 00:11:33 2004 UTC (19 years, 9 months ago) by tim
File size:	15349 byte(s)
Log Message:	adding DUFF parser
#	User	Rev	Content
1	tim	1741	/*
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	gezelter	1490
26	tim	1741	#include "UseTheForce/DUFF.hpp"
27	gezelter	1634	#include "UseTheForce/DarkSide/lj_interface.h"
28	gezelter	1490
29	tim	1741	namespace oopse {
30	gezelter	1490
31	tim	1741	//definition of createDUFF
32			ForceField* createDUFF() {
33			return new DUFF();
34			}
35	gezelter	1490
36	tim	1741	//register createDUFF to ForceFieldCreator
37			ForceFieldFactory::getInstance()->registerForceField("DUFF", createDUFF);
38	gezelter	1490
39	tim	1757
40			ParseState DUFF::getSection(const std::string& section) {
41			ParseState result;
42
43			switch(section) {
44			case "AtomTypes" :
45			result = DUFF::AtomTypeSection;
46			break;
47
48			case "BondTypes" :
49			result = DUFF::BondTypeSection;
50			break;
51
52			case "BendTypes" :
53			result = DUFF::BendTypeSection;
54			break;
55
56			case "TorsionTypes" :
57			result = DUFF::TorsionTypeSection;
58			break;
59			default:
60			result = DUFF::UnknownSection;
61			}
62
63			return result;
64			}
65
66	tim	1741	void DUFF::parse(const std::string& filename) {
67			ifstrstream* ffStream;
68			ffStream = openForceFiledFile(filename);
69			const int bufferSize = 65535;
70	tim	1757	std::string line;
71			char buffer[bufferSize];
72	tim	1741	int lineNo = 0;
73	tim	1757	ParseState currentSection = DUFF::UnknownSection;
74
75			while(ffStream.getline(buffer, bufferSize)){
76	tim	1741	++lineNo;
77	gezelter	1490
78	tim	1757	line = trimSpaces(buffer);
79			//a line begins with "//" is comment
80			if ( line.empty() \|\| (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
81			continue;
82			} else {
83	gezelter	1490
84	tim	1757	switch(currentSection) {
85			case DUFF::AtomTypeSection :
86			parseAtomType(line, lineNo);
87			break;
88	gezelter	1490
89	tim	1757	case DUFF::BondTypeSection :
90			parseBondType(line, lineNo);
91			break;
92	gezelter	1490
93	tim	1757	case DUFF::BendTypeSection :
94			parseBendType(line, lineNo);
95			break;
96
97			case DUFF::TorsionTypeSection :
98			parseTorsionType(line, lineNo);
99			break;
100
101			case DUFF::UnknownSection:
102			StringTokenizer tokenizer(line);
103
104			std::string keyword = tokenizer.nextToken();
105			std::string section = tokenizer.nextToken();
106
107			ParseState newSection = getSection(section);
108			if (keyword != "begin" \|\| keyword != "end") {
109			std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
110			} else if (keyword == "begin") {
111			if (newSection == DUFF::UnknownSection) {
112			std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
113			} else {
114			//enter a new section
115			currentSection = newSection;
116			}
117
118			} else if (keyword == "end"){
119			if (currentSection == newSection) ) {
120			//leave a section
121			currentSection = DUFF::UnknownSection;
122			} else {
123			std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
124			}
125
126			}
127			break;
128			default :
129	tim	1741
130			}
131	gezelter	1653
132	gezelter	1634	}
133	gezelter	1490	}
134
135	tim	1741	delete ffStream;
136	gezelter	1490	}
137
138
139	tim	1757	void DUFF::parseAtomType(const std::string& line, int lineNo){
140			StringTokenizer tokenizer(line);
141
142			if (tokenizer.countTokens() >= 4) {
143			atomTypeName = tokenizer.nextToken();
144			mass = tokenizer.nextTokenAsFloat();
145			epsilon = tokenizer.nextTokenAsFloat();
146			sigma = tokenizer.nextTokenAsFloat();
147
148			atomType = new AtomType();
149			atomType->setName(atomTypeName);
150			atomType->setMass(mass);
151
152			//by default, all of the properties in AtomTypeProperties is set to 0
153			//In Lennard-Jones Force Field, we only need to set Lennard-Jones to true
154			atomType->setLennardJones();
155
156			atomType->setIdent(ident);
157			atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));
158			atomType->addProperty(new DoubleGenericData("Sigma", sigma));
159			atomType->complete();
160			//notify a new LJtype atom type is created
161			newLJtype(&ident, &sigma, &epsilon, &status);
162
163			//add atom type to AtomTypeContainer
164			addAtomType(atomTypeName, atomType);
165			++ident;
166
167			} else {
168			sprintf( painCave.errMsg,
169			"Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
170			"in line %d : %s\n",
171			filename.c_str(), lineNo, line);
172			painCave.severity = OOPSE_ERROR;
173			painCave.isFatal = 1;
174			simError();
175			}
176	gezelter	1490	}
177
178	tim	1757	void DUFF::parseBondType(const std::string& line, int lineNo){
179
180			StringTokenizer tokenizer(line);
181			std::string at1;
182			std::string at2;
183			std::string bt;
184			BondType* bondType = NULL;
185			double b0;
186
187			int nTokens = tokenizer.countTokens();
188
189			if (nTokens < 4) {
190
191			return;
192			}
193
194			at1 = tokenizer.nextToken();
195			at2 = tokenizer.nextToken();
196			bt = tokenizer.nextToken();
197			b0 = tokenizer.nextTokenAsDouble();
198			nTokens -= 4;
199
200			//switch is a maintain nightmare
201			switch(bt) {
202			case "FixedBondType" :
203			bondType = new FixedBondType();
204			break;
205
206			case "HarmonicBondType" :
207			if (nTokens < 1) {
208
209			} else {
210
211			double kb = tokenizer.nextTokenAsDouble();
212			bondType = new HarmonicBondType(b0, kb);
213			}
214
215			break;
216
217			case "CubicBondType" :
218			if (nTokens < 4) {
219
220			} else {
221
222			double k3 = tokenizer.nextTokenAsDouble();
223			double k2 = tokenizer.nextTokenAsDouble();
224			double k1 = tokenizer.nextTokenAsDouble();
225			double k0 = tokenizer.nextTokenAsDouble();
226
227			bondType = new CubicBondType(b0, k3, k2, k1, k0);
228			}
229			break;
230
231			case "QuadraticBondType" :
232			if (nTokens < 5) {
233
234			} else {
235
236			b0 = tokenizer.nextTokenAsDouble();
237			double k4 = tokenizer.nextTokenAsDouble();
238			double k3 = tokenizer.nextTokenAsDouble();
239			double k2 = tokenizer.nextTokenAsDouble();
240			double k1 = tokenizer.nextTokenAsDouble();
241			double k0 = tokenizer.nextTokenAsDouble();
242
243			bondType = new QuadraticBondType(b0, k4, k3, k2, k1, k0);
244			}
245			break;
246
247			case "PolynomialBondType " :
248			if (nTokens < 2 \|\| nTokens % 2 != 0) {
249
250			} else {
251			int nPairs = nTokens / 2;
252			int power;
253			double coefficient;
254			PolynomialBondType pbt = new PolynomialBondType();
255
256			for (int i = 0; i < nPairs; ++i) {
257			power = tokenizer.nextTokenAsInt();
258			coefficient = tokenizer.nextTokenAsDouble();
259
260			if (power < 0) {
261
262			} else {
263			pbt->setCoefficient(power, coefficient);
264			}
265			}
266			}
267
268			break;
269
270			default:
271
272			}
273
274			if (bondType != NULL) {
275			addBondType(at1, at2, bondType);
276			}
277	gezelter	1490	}
278
279	tim	1757	void DUFF::parseBendType(const std::string& line, int lineNo){
280			StringTokenizer tokenizer(line);
281			std::string at1;
282			std::string at2;
283			std::string at3;
284			std::string bt;
285			double theta0;
286			BendType* bendType = NULL;
287
288			int nTokens = tokenizer.countTokens();
289
290			if (nTokens < 5) {
291
292			return;
293			}
294
295			at1 = tokenizer.nextToken();
296			at2 = tokenizer.nextToken();
297			at3 = tokenizer.nextToken();
298			bt = tokenizer.nextToken();
299			theta0 = tokenizer.nextTokenAsDouble();
300			nTokens -= 5;
301
302			//switch is a maintain nightmare
303			switch(bt) {
304
305			case "HarmonicBendType" :
306
307			if (nTokens < 1) {
308
309			} else {
310
311			double ktheta = tokenizer.nextTokenAsDouble();
312			bendType = new HarmonicBendType(theta0, ktheta);
313			}
314			break;
315			case "GhostBendType" :
316			if (nTokens < 1) {
317
318			} else {
319			double ktheta = tokenizer.nextTokenAsDouble();
320			bendType = new HarmonicBendType(theta0, ktheta);
321			}
322			break;
323
324			case "UreyBradleyBendType" :
325			if (nTokens < 3) {
326
327			} else {
328			double ktheta = tokenizer.nextTokenAsDouble();
329			double s0 = tokenizer.nextTokenAsDouble();
330			double kub = tokenizer.nextTokenAsDouble();
331			bendType = new UreyBradleyBendType(theta0, ktheta, s0, kub);
332			}
333			break;
334
335			case "CubicBendType" :
336			if (nTokens < 4) {
337
338			} else {
339
340			double k3 = tokenizer.nextTokenAsDouble();
341			double k2 = tokenizer.nextTokenAsDouble();
342			double k1 = tokenizer.nextTokenAsDouble();
343			double k0 = tokenizer.nextTokenAsDouble();
344
345			bendType = new CubicBendType(theta0, k3, k2, k1, k0);
346			}
347			break;
348
349			case "QuadraticBendType" :
350			if (nTokens < 5) {
351
352			} else {
353
354			theta0 = tokenizer.nextTokenAsDouble();
355			double k4 = tokenizer.nextTokenAsDouble();
356			double k3 = tokenizer.nextTokenAsDouble();
357			double k2 = tokenizer.nextTokenAsDouble();
358			double k1 = tokenizer.nextTokenAsDouble();
359			double k0 = tokenizer.nextTokenAsDouble();
360
361			bendType = new QuadraticBendType(theta0, k4, k3, k2, k1, k0);
362			}
363			break;
364
365			case "PolynomialBendType " :
366			if (nTokens < 2 \|\| nTokens % 2 != 0) {
367
368			} else {
369			int nPairs = nTokens / 2;
370			int power;
371			double coefficient;
372			PolynomialBendType* pbt = new PolynomialBendType();
373
374			for (int i = 0; i < nPairs; ++i) {
375			power = tokenizer.nextTokenAsInt();
376			coefficient = tokenizer.nextTokenAsDouble();
377
378			if (power < 0) {
379
380			} else {
381			pbt->setCoefficient(power, coefficient);
382			}
383			}
384			}
385
386			break;
387
388			default:
389
390			}
391
392			if (bendType != NULL) {
393			addBendType(at1, at2, at3, bendType);
394			}
395
396	gezelter	1490	}
397
398	tim	1757	void DUFF::parseTorsionType(const std::string& line, int lineNo){
399			StringTokenizer tokenizer(line);
400			std::string at1;
401			std::string at2;
402			std::string at3;
403			std::string at4;
404			std::string tt;
405			TorsionType* bendType = NULL;
406	gezelter	1490
407	tim	1757	int nTokens = tokenizer.countTokens();
408
409			if (nTokens < 5) {
410
411			return;
412			}
413
414			at1 = tokenizer.nextToken();
415			at2 = tokenizer.nextToken();
416			at3 = tokenizer.nextToken();
417			at4 = tokenizer.nextToken();
418			tt = tokenizer.nextToken();
419
420			nTokens -= 5;
421
422			switch(tt) {
423
424			case "CubicTorsionType" :
425			if (nTokens < 4) {
426
427			} else {
428
429			double k3 = tokenizer.nextTokenAsDouble();
430			double k2 = tokenizer.nextTokenAsDouble();
431			double k1 = tokenizer.nextTokenAsDouble();
432			double k0 = tokenizer.nextTokenAsDouble();
433
434			bendType = new CubicTorsionType(k3, k2, k1, k0);
435			}
436			break;
437
438			case "QuadraticTorsionType" :
439			if (nTokens < 5) {
440
441			} else {
442
443			theta0 = tokenizer.nextTokenAsDouble();
444			double k4 = tokenizer.nextTokenAsDouble();
445			double k3 = tokenizer.nextTokenAsDouble();
446			double k2 = tokenizer.nextTokenAsDouble();
447			double k1 = tokenizer.nextTokenAsDouble();
448			double k0 = tokenizer.nextTokenAsDouble();
449
450			bendType = new QuadraticTorsionType( k4, k3, k2, k1, k0);
451			}
452			break;
453
454			case "PolynomialTorsionType " :
455			if (nTokens < 2 \|\| nTokens % 2 != 0) {
456
457			} else {
458			int nPairs = nTokens / 2;
459			int power;
460			double coefficient;
461			PolynomialTorsionType* pbt = new PolynomialTorsionType();
462
463			for (int i = 0; i < nPairs; ++i) {
464			power = tokenizer.nextTokenAsInt();
465			coefficient = tokenizer.nextTokenAsDouble();
466
467			if (power < 0) {
468
469			} else {
470			pbt->setCoefficient(power, coefficient);
471			}
472			}
473			}
474
475			break;
476			case "CharmmTorsionType" :
477
478			if (nTokens < 3 \|\| nTokens % 3 != 0) {
479
480			} else {
481			int nSets = nTokens / 3;
482
483			double kchi;
484			int n;
485			double delta;
486			CharmmTorsionType* ctt = new CharmmTorsionType();
487
488			for (int i = 0; i < nSets; ++i) {
489			kchi = tokenizer.nextTokenAsDouble();
490			n = tokenizer.nextTokenAsInt();
491			delta = tokenizer.nextTokenAsDouble();
492
493			if (n < 0 \|\| n > 6) {
494
495			} else {
496			ctt->setCharmmTorsionParameter(kchi, n, delta);
497			}
498			}
499			}
500			default:
501
502			}
503
504			if (bendType != NULL) {
505			addTorsionType(at1, at2, at3, bendType);
506			}
507	gezelter	1490	}
508
509	tim	1741	} //end namespace oopse
510	gezelter	1490