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"
#include "UseTheForce/DarkSide/charge_interface.h"
#include "UseTheForce/DarkSide/dipole_interface.h"
#include "UseTheForce/DarkSide/sticky_interface.h"

namespace oopse {

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

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


ParseState DUFF::getSection(const std::string& section) {
    ParseState result;
    
    switch(section) {
        case "AtomTypes" :
            result = DUFF::AtomTypeSection;
            break;
        case "DirectionalAtomTypes" :
            result = DUFF::DirectionalAtomTypeSection;
            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;
    int atomIdent = getNAtomType() + 1;  //atom's indent begins from 1 (since only fortran's array begins from 1)
    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, atomIdent);
                    break;

                case DUFF::DirectionalAtomTypeSection :
                    parseDirectionalAtomType(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, int& ident){ 
    StringTokenizer tokenizer(line);
    int nTokens = tokenizer.countTokens();    

    //in AtomTypeSection, a line at least contains 5 tokens
    //atomTypeName, is Directional, isLJ, isCharge and mass
    if (nTokens < 5)  {
                      
    } else {

        std::string atomTypeName = tokenizer.nextToken();
        bool isDirectional = tokenizer.nextTokenAsBool();
        bool isLJ = tokenizer.nextTokenAsBool();
        bool isCharge = tokenizer.nextTokenAsBool();        
        double mass = tokenizer.nextTokenAsDouble();
        double epsilon;
        double sigma;
        double charge;
        nTokens -= 5;

        //parse epsilon and sigma
        if (isLJ) {
            if (nTokens >= 2) {
                epsilon = tokenizer.nextTokenAsDouble();
                sigma = tokenizer.nextTokenAsDouble();
                nTokens -= 2;
            } else {

            }
        }

        //parse charge 
        if (isCharge) {
            if (nTokens >= 1) {
                charge = tokenizer.nextTokenAsDouble();
                nTokens -= 1;
            } else {

            }
        }
        
        AtomType* atomType;
        if (isDirectional) {
            atomType = new DirectionalAtomType();        
        } else {
            atomType = new AtomType();
        }
        
        atomType->setName(atomTypeName);
        atomType->setMass(mass);
        
        if (isLJ) {
            atomType->setLennardJones();
        }

        if (isCharge) {
            atomType->setCharge();
        }
        
        atomType->setIdent(ident); 
        
        atomType->complete();

        int setLJStatus;
        
        //notify a new LJtype atom type is created
        if (isLJ) {
            newLJtype(&ident, &sigma, &epsilon, &setLJStatus);
        }

        int setChargeStatus;
        if (isCharge) {
            newChargeType(&ident, &charge, &setChargeStatus)
        }

        if (setLJStatus && setChargeStatus) {
            //add atom type to AtomTypeContainer
            addAtomType(atomTypeName, atomType);
            ++ident;
        } else {
            //error in notifying fortran
            delete atomType;
        }
    }    
    
}


void DUFF::parseDirectionalAtomType(const std::string& line, int lineNo) {
    StringTokenizer tokenizer(line);
    int nTokens = tokenizer.countTokens();  
    
    //in DirectionalAtomTypeSection, a line at least contains 6 tokens
    //AtomTypeName, isDipole, isSticky, I_xx, I_yy and I_zz
    if (nTokens < 6) {
        std::cerr << "Not enought tokens" << std::endl;
    } else {


        std::string atomTypeName = tokenizer.nextToken();
        bool isDipole = tokenizer.nextTokenAsBool();
        bool isSticky = tokenizer.nextTokenAsBool();
        double Ixx = tokenizer.nextTokenAsDouble();
        double Iyy = tokenizer.nextTokenAsDouble();
        double Izz = tokenizer.nextTokenAsDouble();
        nTokens -= 6;

        AtomType* atomType = getAtomType(atomTypeName);
        if (atomType == NULL) {

        }
        
        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
        if (dAtomType == NULL) {


        }

        if (isDipole) {
            dAtomType->setDipole();
        }

        if (isSticky) {
            dAtomType->setSticky();
        }
        
        Mat3x3d inertialMat;
        inertialMat(0, 0) = Ixx;
        inertialMat(1, 1) = Iyy;
        inertialMat(2, 2) = Izz;        
        dAtomType->setI(inertialMat);

        //read dipole moment
        double dipole;
        if (isDipole) {
            if (nTokens >= 1) {
                dipole = tokenizer.nextTokenAsDouble();
                nTokens -= 1;
            } else {

            }
        }

        //read sticky parameters
        double w0;
        double v0;
        double v0p;
        double rl;
        double ru;
        double rlp;
        double rup;        
        if (isSticky) {
            if (nTokens >= 7) {
                w0 = tokenizer.nextTokenAsDouble();
                v0 = tokenizer.nextTokenAsDouble();
                v0p = tokenizer.nextTokenAsDouble();
                rl = tokenizer.nextTokenAsDouble();
                ru = tokenizer.nextTokenAsDouble();
                rlp = tokenizer.nextTokenAsDouble();
                rup = tokenizer.nextTokenAsDouble();     
                nTokens -= 7;
            } else {

            }
        }


        //notify fotran a newDipoleType is created
        int ident = dAtomType->getIdent();
        int setDipoleStatus;
        if (isDipole) {
            newDipoleType(&ident, &dipole, &setDipoleStatus);
        }               

        //notify fotran a StickyType is created
        int setStickyStatus;
        if (isSticky) {
            makeStickyType( &w0, &v0, &v0p, &rl, &ru, &rlp, &rup);
        }


        if (!setDipoleStatus || !setStickyStatus) {

        }
        
    }
}

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 "Fixed" :
            bondType = new FixedBondType();
            break;
            
        case "Harmonic" :
            if (nTokens < 1) {

            } else {

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

            break;

        case "Cubic" :
            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 "Quartic" :
            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 "Polynomial" :
            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();
                    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 "Harmonic" :
            
            if (nTokens < 1) {

            } else {

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

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

        case "UreyBradley" :
            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 "Cubic" :
            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 "Quartic" :
            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 "Polynomial" :
            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();
                    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* torsionType = 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 "Cubic" :
            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 "Quartic" :
            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 "Polynomial" :
            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();
                    pbt->setCoefficient(power, coefficient);
                }
            }
            
            break;
        case "Charmm" :
            
            if (nTokens < 3 || nTokens % 3 != 0) {

            } else {
                int nSets = nTokens / 3;
  
                CharmmTorsionType* ctt = new CharmmTorsionType();
                
                for (int i = 0; i < nSets; ++i) {
                    double kchi = tokenizer.nextTokenAsDouble();
                    int n = tokenizer.nextTokenAsInt();
                    double delta = tokenizer.nextTokenAsDouble();
    
                    ctt->setCharmmTorsionParameter(kchi, n, delta);
                }
            }
        default:
            
    }

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

} //end namespace oopse
Revision:	1767
Committed:	Mon Nov 22 23:04:16 2004 UTC (19 years, 9 months ago) by tim
File size:	18851 byte(s)
Log Message:	minor change in DUFF and EAM
#	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	tim	1758	#include "UseTheForce/DarkSide/charge_interface.h"
29			#include "UseTheForce/DarkSide/dipole_interface.h"
30			#include "UseTheForce/DarkSide/sticky_interface.h"
31	gezelter	1490
32	tim	1741	namespace oopse {
33	gezelter	1490
34	tim	1741	//definition of createDUFF
35			ForceField* createDUFF() {
36			return new DUFF();
37			}
38	gezelter	1490
39	tim	1758	//register createDUFF to ForceFieldFactory
40	tim	1741	ForceFieldFactory::getInstance()->registerForceField("DUFF", createDUFF);
41	gezelter	1490
42	tim	1757
43			ParseState DUFF::getSection(const std::string& section) {
44			ParseState result;
45
46			switch(section) {
47			case "AtomTypes" :
48			result = DUFF::AtomTypeSection;
49			break;
50	tim	1758	case "DirectionalAtomTypes" :
51			result = DUFF::DirectionalAtomTypeSection;
52			break;
53	tim	1757
54			case "BondTypes" :
55			result = DUFF::BondTypeSection;
56			break;
57
58			case "BendTypes" :
59			result = DUFF::BendTypeSection;
60			break;
61
62			case "TorsionTypes" :
63			result = DUFF::TorsionTypeSection;
64			break;
65			default:
66			result = DUFF::UnknownSection;
67			}
68
69			return result;
70			}
71
72	tim	1741	void DUFF::parse(const std::string& filename) {
73			ifstrstream* ffStream;
74			ffStream = openForceFiledFile(filename);
75			const int bufferSize = 65535;
76	tim	1757	std::string line;
77			char buffer[bufferSize];
78	tim	1741	int lineNo = 0;
79	tim	1760	int atomIdent = getNAtomType() + 1; //atom's indent begins from 1 (since only fortran's array begins from 1)
80	tim	1757	ParseState currentSection = DUFF::UnknownSection;
81
82			while(ffStream.getline(buffer, bufferSize)){
83	tim	1741	++lineNo;
84	gezelter	1490
85	tim	1757	line = trimSpaces(buffer);
86			//a line begins with "//" is comment
87			if ( line.empty() \|\| (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
88			continue;
89			} else {
90	gezelter	1490
91	tim	1757	switch(currentSection) {
92			case DUFF::AtomTypeSection :
93	tim	1759	parseAtomType(line, lineNo, atomIdent);
94	tim	1757	break;
95	gezelter	1490
96	tim	1759	case DUFF::DirectionalAtomTypeSection :
97			parseDirectionalAtomType(line, lineNo);
98			break;
99
100	tim	1757	case DUFF::BondTypeSection :
101			parseBondType(line, lineNo);
102			break;
103	gezelter	1490
104	tim	1757	case DUFF::BendTypeSection :
105			parseBendType(line, lineNo);
106			break;
107
108			case DUFF::TorsionTypeSection :
109			parseTorsionType(line, lineNo);
110			break;
111
112			case DUFF::UnknownSection:
113			StringTokenizer tokenizer(line);
114
115			std::string keyword = tokenizer.nextToken();
116			std::string section = tokenizer.nextToken();
117
118			ParseState newSection = getSection(section);
119			if (keyword != "begin" \|\| keyword != "end") {
120			std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
121			} else if (keyword == "begin") {
122			if (newSection == DUFF::UnknownSection) {
123			std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
124			} else {
125			//enter a new section
126			currentSection = newSection;
127			}
128
129			} else if (keyword == "end"){
130			if (currentSection == newSection) ) {
131			//leave a section
132			currentSection = DUFF::UnknownSection;
133			} else {
134			std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
135			}
136
137			}
138			break;
139			default :
140	tim	1741
141			}
142	gezelter	1653
143	gezelter	1634	}
144	gezelter	1490	}
145
146	tim	1741	delete ffStream;
147	gezelter	1490	}
148
149	tim	1758	void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){
150	tim	1757	StringTokenizer tokenizer(line);
151	tim	1758	int nTokens = tokenizer.countTokens();
152	tim	1757
153	tim	1758	//in AtomTypeSection, a line at least contains 5 tokens
154			//atomTypeName, is Directional, isLJ, isCharge and mass
155			if (nTokens < 5) {
156	tim	1759
157	tim	1758	} else {
158
159			std::string atomTypeName = tokenizer.nextToken();
160			bool isDirectional = tokenizer.nextTokenAsBool();
161			bool isLJ = tokenizer.nextTokenAsBool();
162			bool isCharge = tokenizer.nextTokenAsBool();
163			double mass = tokenizer.nextTokenAsDouble();
164			double epsilon;
165			double sigma;
166			double charge;
167			nTokens -= 5;
168
169			//parse epsilon and sigma
170			if (isLJ) {
171			if (nTokens >= 2) {
172			epsilon = tokenizer.nextTokenAsDouble();
173			sigma = tokenizer.nextTokenAsDouble();
174			nTokens -= 2;
175			} else {
176
177			}
178			}
179
180			//parse charge
181			if (isCharge) {
182			if (nTokens >= 1) {
183			charge = tokenizer.nextTokenAsDouble();
184			nTokens -= 1;
185			} else {
186
187			}
188			}
189
190			AtomType* atomType;
191			if (isDirectional) {
192			atomType = new DirectionalAtomType();
193			} else {
194			atomType = new AtomType();
195			}
196
197	tim	1757	atomType->setName(atomTypeName);
198			atomType->setMass(mass);
199
200	tim	1760	if (isLJ) {
201			atomType->setLennardJones();
202			}
203	tim	1757
204	tim	1760	if (isCharge) {
205			atomType->setCharge();
206			}
207
208	tim	1758	atomType->setIdent(ident);
209
210	tim	1757	atomType->complete();
211	tim	1758
212			int setLJStatus;
213
214	tim	1757	//notify a new LJtype atom type is created
215	tim	1758	if (isLJ) {
216			newLJtype(&ident, &sigma, &epsilon, &setLJStatus);
217			}
218	tim	1757
219	tim	1758	int setChargeStatus;
220			if (isCharge) {
221			newChargeType(&ident, &charge, &setChargeStatus)
222			}
223
224			if (setLJStatus && setChargeStatus) {
225			//add atom type to AtomTypeContainer
226			addAtomType(atomTypeName, atomType);
227			++ident;
228			} else {
229			//error in notifying fortran
230			delete atomType;
231			}
232	tim	1757	}
233	tim	1758
234	gezelter	1490	}
235
236	tim	1758
237			void DUFF::parseDirectionalAtomType(const std::string& line, int lineNo) {
238			StringTokenizer tokenizer(line);
239			int nTokens = tokenizer.countTokens();
240
241	tim	1765	//in DirectionalAtomTypeSection, a line at least contains 6 tokens
242	tim	1758	//AtomTypeName, isDipole, isSticky, I_xx, I_yy and I_zz
243			if (nTokens < 6) {
244	tim	1765	std::cerr << "Not enought tokens" << std::endl;
245	tim	1758	} else {
246
247
248			std::string atomTypeName = tokenizer.nextToken();
249			bool isDipole = tokenizer.nextTokenAsBool();
250			bool isSticky = tokenizer.nextTokenAsBool();
251			double Ixx = tokenizer.nextTokenAsDouble();
252			double Iyy = tokenizer.nextTokenAsDouble();
253			double Izz = tokenizer.nextTokenAsDouble();
254			nTokens -= 6;
255
256			AtomType* atomType = getAtomType(atomTypeName);
257			if (atomType == NULL) {
258
259			}
260
261			DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
262			if (dAtomType == NULL) {
263
264
265			}
266
267	tim	1760	if (isDipole) {
268			dAtomType->setDipole();
269			}
270	tim	1758
271	tim	1760	if (isSticky) {
272			dAtomType->setSticky();
273			}
274
275	tim	1758	Mat3x3d inertialMat;
276			inertialMat(0, 0) = Ixx;
277	tim	1765	inertialMat(1, 1) = Iyy;
278			inertialMat(2, 2) = Izz;
279	tim	1758	dAtomType->setI(inertialMat);
280
281			//read dipole moment
282			double dipole;
283			if (isDipole) {
284			if (nTokens >= 1) {
285			dipole = tokenizer.nextTokenAsDouble();
286			nTokens -= 1;
287			} else {
288
289			}
290			}
291
292			//read sticky parameters
293			double w0;
294			double v0;
295			double v0p;
296			double rl;
297			double ru;
298			double rlp;
299			double rup;
300			if (isSticky) {
301			if (nTokens >= 7) {
302			w0 = tokenizer.nextTokenAsDouble();
303			v0 = tokenizer.nextTokenAsDouble();
304			v0p = tokenizer.nextTokenAsDouble();
305			rl = tokenizer.nextTokenAsDouble();
306			ru = tokenizer.nextTokenAsDouble();
307			rlp = tokenizer.nextTokenAsDouble();
308			rup = tokenizer.nextTokenAsDouble();
309			nTokens -= 7;
310			} else {
311
312			}
313			}
314
315
316			//notify fotran a newDipoleType is created
317			int ident = dAtomType->getIdent();
318			int setDipoleStatus;
319			if (isDipole) {
320			newDipoleType(&ident, &dipole, &setDipoleStatus);
321			}
322
323			//notify fotran a StickyType is created
324			int setStickyStatus;
325			if (isSticky) {
326			makeStickyType( &w0, &v0, &v0p, &rl, &ru, &rlp, &rup);
327			}
328
329
330			if (!setDipoleStatus \|\| !setStickyStatus) {
331
332			}
333
334			}
335			}
336
337	tim	1757	void DUFF::parseBondType(const std::string& line, int lineNo){
338
339			StringTokenizer tokenizer(line);
340			std::string at1;
341			std::string at2;
342			std::string bt;
343			BondType* bondType = NULL;
344			double b0;
345
346			int nTokens = tokenizer.countTokens();
347
348			if (nTokens < 4) {
349
350			return;
351			}
352
353			at1 = tokenizer.nextToken();
354			at2 = tokenizer.nextToken();
355			bt = tokenizer.nextToken();
356			b0 = tokenizer.nextTokenAsDouble();
357			nTokens -= 4;
358
359			//switch is a maintain nightmare
360			switch(bt) {
361	tim	1767	case "Fixed" :
362	tim	1757	bondType = new FixedBondType();
363			break;
364
365	tim	1767	case "Harmonic" :
366	tim	1757	if (nTokens < 1) {
367
368			} else {
369
370			double kb = tokenizer.nextTokenAsDouble();
371			bondType = new HarmonicBondType(b0, kb);
372			}
373
374			break;
375
376	tim	1767	case "Cubic" :
377	tim	1757	if (nTokens < 4) {
378
379			} else {
380
381			double k3 = tokenizer.nextTokenAsDouble();
382			double k2 = tokenizer.nextTokenAsDouble();
383			double k1 = tokenizer.nextTokenAsDouble();
384			double k0 = tokenizer.nextTokenAsDouble();
385
386			bondType = new CubicBondType(b0, k3, k2, k1, k0);
387			}
388			break;
389
390	tim	1767	case "Quartic" :
391	tim	1757	if (nTokens < 5) {
392
393			} else {
394
395			b0 = tokenizer.nextTokenAsDouble();
396			double k4 = tokenizer.nextTokenAsDouble();
397			double k3 = tokenizer.nextTokenAsDouble();
398			double k2 = tokenizer.nextTokenAsDouble();
399			double k1 = tokenizer.nextTokenAsDouble();
400			double k0 = tokenizer.nextTokenAsDouble();
401
402			bondType = new QuadraticBondType(b0, k4, k3, k2, k1, k0);
403			}
404			break;
405
406	tim	1767	case "Polynomial" :
407	tim	1757	if (nTokens < 2 \|\| nTokens % 2 != 0) {
408
409			} else {
410			int nPairs = nTokens / 2;
411			int power;
412			double coefficient;
413			PolynomialBondType pbt = new PolynomialBondType();
414
415			for (int i = 0; i < nPairs; ++i) {
416			power = tokenizer.nextTokenAsInt();
417			coefficient = tokenizer.nextTokenAsDouble();
418	tim	1758	pbt->setCoefficient(power, coefficient);
419	tim	1757	}
420			}
421
422			break;
423
424			default:
425
426			}
427
428			if (bondType != NULL) {
429			addBondType(at1, at2, bondType);
430			}
431	gezelter	1490	}
432
433	tim	1757	void DUFF::parseBendType(const std::string& line, int lineNo){
434			StringTokenizer tokenizer(line);
435			std::string at1;
436			std::string at2;
437			std::string at3;
438			std::string bt;
439			double theta0;
440			BendType* bendType = NULL;
441
442			int nTokens = tokenizer.countTokens();
443
444			if (nTokens < 5) {
445
446			return;
447			}
448
449			at1 = tokenizer.nextToken();
450			at2 = tokenizer.nextToken();
451			at3 = tokenizer.nextToken();
452			bt = tokenizer.nextToken();
453			theta0 = tokenizer.nextTokenAsDouble();
454			nTokens -= 5;
455
456			//switch is a maintain nightmare
457			switch(bt) {
458
459	tim	1767	case "Harmonic" :
460	tim	1757
461			if (nTokens < 1) {
462
463			} else {
464
465			double ktheta = tokenizer.nextTokenAsDouble();
466			bendType = new HarmonicBendType(theta0, ktheta);
467			}
468			break;
469	tim	1767	case "GhostBend" :
470	tim	1757	if (nTokens < 1) {
471
472			} else {
473			double ktheta = tokenizer.nextTokenAsDouble();
474			bendType = new HarmonicBendType(theta0, ktheta);
475			}
476			break;
477
478	tim	1767	case "UreyBradley" :
479	tim	1757	if (nTokens < 3) {
480
481			} else {
482			double ktheta = tokenizer.nextTokenAsDouble();
483			double s0 = tokenizer.nextTokenAsDouble();
484			double kub = tokenizer.nextTokenAsDouble();
485			bendType = new UreyBradleyBendType(theta0, ktheta, s0, kub);
486			}
487			break;
488
489	tim	1767	case "Cubic" :
490	tim	1757	if (nTokens < 4) {
491
492			} else {
493
494			double k3 = tokenizer.nextTokenAsDouble();
495			double k2 = tokenizer.nextTokenAsDouble();
496			double k1 = tokenizer.nextTokenAsDouble();
497			double k0 = tokenizer.nextTokenAsDouble();
498
499			bendType = new CubicBendType(theta0, k3, k2, k1, k0);
500			}
501			break;
502
503	tim	1767	case "Quartic" :
504	tim	1757	if (nTokens < 5) {
505
506			} else {
507
508			theta0 = tokenizer.nextTokenAsDouble();
509			double k4 = tokenizer.nextTokenAsDouble();
510			double k3 = tokenizer.nextTokenAsDouble();
511			double k2 = tokenizer.nextTokenAsDouble();
512			double k1 = tokenizer.nextTokenAsDouble();
513			double k0 = tokenizer.nextTokenAsDouble();
514
515			bendType = new QuadraticBendType(theta0, k4, k3, k2, k1, k0);
516			}
517			break;
518
519	tim	1767	case "Polynomial" :
520	tim	1757	if (nTokens < 2 \|\| nTokens % 2 != 0) {
521
522			} else {
523			int nPairs = nTokens / 2;
524			int power;
525			double coefficient;
526			PolynomialBendType* pbt = new PolynomialBendType();
527
528			for (int i = 0; i < nPairs; ++i) {
529			power = tokenizer.nextTokenAsInt();
530			coefficient = tokenizer.nextTokenAsDouble();
531	tim	1758	pbt->setCoefficient(power, coefficient);
532	tim	1757	}
533			}
534
535			break;
536
537			default:
538
539			}
540
541			if (bendType != NULL) {
542			addBendType(at1, at2, at3, bendType);
543			}
544
545	gezelter	1490	}
546
547	tim	1757	void DUFF::parseTorsionType(const std::string& line, int lineNo){
548			StringTokenizer tokenizer(line);
549			std::string at1;
550			std::string at2;
551			std::string at3;
552			std::string at4;
553			std::string tt;
554	tim	1767	TorsionType* torsionType = NULL;
555	gezelter	1490
556	tim	1757	int nTokens = tokenizer.countTokens();
557
558			if (nTokens < 5) {
559
560			return;
561			}
562
563			at1 = tokenizer.nextToken();
564			at2 = tokenizer.nextToken();
565			at3 = tokenizer.nextToken();
566			at4 = tokenizer.nextToken();
567			tt = tokenizer.nextToken();
568
569			nTokens -= 5;
570
571			switch(tt) {
572
573	tim	1767	case "Cubic" :
574	tim	1757	if (nTokens < 4) {
575
576			} else {
577
578			double k3 = tokenizer.nextTokenAsDouble();
579			double k2 = tokenizer.nextTokenAsDouble();
580			double k1 = tokenizer.nextTokenAsDouble();
581			double k0 = tokenizer.nextTokenAsDouble();
582
583			bendType = new CubicTorsionType(k3, k2, k1, k0);
584			}
585			break;
586
587	tim	1767	case "Quartic" :
588	tim	1757	if (nTokens < 5) {
589
590			} else {
591
592			theta0 = tokenizer.nextTokenAsDouble();
593			double k4 = tokenizer.nextTokenAsDouble();
594			double k3 = tokenizer.nextTokenAsDouble();
595			double k2 = tokenizer.nextTokenAsDouble();
596			double k1 = tokenizer.nextTokenAsDouble();
597			double k0 = tokenizer.nextTokenAsDouble();
598
599			bendType = new QuadraticTorsionType( k4, k3, k2, k1, k0);
600			}
601			break;
602
603	tim	1767	case "Polynomial" :
604	tim	1757	if (nTokens < 2 \|\| nTokens % 2 != 0) {
605
606			} else {
607			int nPairs = nTokens / 2;
608			int power;
609			double coefficient;
610			PolynomialTorsionType* pbt = new PolynomialTorsionType();
611
612			for (int i = 0; i < nPairs; ++i) {
613			power = tokenizer.nextTokenAsInt();
614			coefficient = tokenizer.nextTokenAsDouble();
615	tim	1758	pbt->setCoefficient(power, coefficient);
616	tim	1757	}
617			}
618
619			break;
620	tim	1767	case "Charmm" :
621	tim	1757
622			if (nTokens < 3 \|\| nTokens % 3 != 0) {
623
624			} else {
625			int nSets = nTokens / 3;
626
627			CharmmTorsionType* ctt = new CharmmTorsionType();
628
629			for (int i = 0; i < nSets; ++i) {
630	tim	1758	double kchi = tokenizer.nextTokenAsDouble();
631			int n = tokenizer.nextTokenAsInt();
632			double delta = tokenizer.nextTokenAsDouble();
633	tim	1757
634	tim	1758	ctt->setCharmmTorsionParameter(kchi, n, delta);
635	tim	1757	}
636			}
637			default:
638
639			}
640
641	tim	1767	if (torsionType != NULL) {
642			addTorsionType(at1, at2, at3, at4, torsionType);
643	tim	1757	}
644	gezelter	1490	}
645
646	tim	1741	} //end namespace oopse